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Abstract 

Two  new  linear  reconstruction  techniques  are  developed  to  improve  the  resolution 
of  images  collected  by  ground-based  telescopes  imaging  through  atmospheric  turbulence. 
The  classical  approach  involves  the  application  of  constrained  least  squares  (CLS)  to  the 
deconvolution  from  wavefront  sensing  (DWFS)  technique.  The  new  algorithm  incorporates 
blur  and  noise  models  to  select  the  appropriate  regularization  constant  automatically.  In  all 
cases  examined,  the  Newton- Raphson  minimization  converged  to  a  solution  in  less  than  10 
iterations.  The  non-iterative  Bayesian  approach  involves  the  development  of  a  new  vector 
Wiener  filter  which  is  optimal  with  respect  to  mean  square  error  (MSE)  for  a  non-stationary 
object  class  degraded  by  atmospheric  turbulence  and  measurement  noise.  This  research 
involves  the  first  extension  of  the  Wiener  filter  to  account  properly  for  shot  noise  and  an 
unknown,  random  optical  transfer  function  (OTF).  The  vector  Wiener  filter  provides  superior 
reconstructions  when  compared  to  the  traditional  scalar  Wiener  filter  for  a  non-stationary 
object  class.  In  addition,  the  new  filter  can  provide  a  superresolution  capability  when  the 
object’s  Fourier  domain  statistics  are  known  for  spatial  frequencies  beyond  the  OTF  cutoff. 
A  generalized  performance  and  robustness  study  of  the  vector  Wiener  filter  showed  that  MSE 
performance  is  fundamentally  limited  by  object  signal-to-noise  ratio  (SNR)  and  correlation 
between  object  pixels. 
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Linear  Reconstruction 


of  Non- Stationary  Image  Ensembles 
Incorporating  Blur  and  Noise  Models 


1.  Introduction 

1.1  Problem  Statement 

Since  before  the  time  of  Galileo  and  Newton,  man  has  used  optical  devices  to  form 
images  of  distant  objects.  The  term  image  refers  to  the  two  dimensional  picture  associated 
with  the  light  or  irradiance  collected  by  an  imaging  system  such  as  the  eye,  a  camera,  or 
telescope.  The  term  object  denotes  the  light  or  radiant  exitance  that  caused  the  image  to  be 
formed.  A  perfect  optical  system  will  produce  an  image  that  is  identical  to  the  object  within 
the  limits  of  diffraction.  In  reality,  a  loss  of  resolution  may  occur  due  to  blur  associated  with 
distortions  in  the  optical  device  or  randomness  in  the  imaging  medium.  An  image  may  also 
be  distorted  by  noise  due  to  low  light  level  or  limitations  in  the  recording  device.  The  imaging 
scenario  describes  the  degradations  affecting  an  optical  system  in  a  given  application.  For 
example,  images  collected  using  a  ground-based  astronomical  telescope  are  degraded  by  the 
turbulent  atmosphere  and  film-grain  or  electronic  detector  noise. 

With  the  widespread  availability  of  computers,  the  concept  of  a  digital  image  has 
become  important.  A  digital  image  is  an  array  of  real  or  complex  numbers  which  represents 
a  sampled  version  of  the  two  dimensional  continuous  image  described  above.  The  elements 
of  an  image  array  are  known  as  pixels.  In  this  mathematical  form,  a  distorted  image  can 
be  manipulated  by  a  computer  using  a  variety  of  techniques  [41].  Image  reconstruction 
refers  to  digital  image  processing  techniques  that  attempt  to  recover  an  accurate  object 
estimate  based  on  a  priori  knowledge  of  the  imaging  scenario.  Statistical  estimation  theory 
plays  an  important  role  in  many  modern  reconstruction  algorithms.  This  theory  can  be 
divided  into  two  main  approaches;  classical  and  Bayesian  [43].  In  a  classical  approach,  the 
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parameter  of  interest  is  assumed  to  be  a  deterministic  but  unknown  constant.  In  contrast, 
the  unknown  parameter  is  assumed  to  be  a  random  variable  in  a  Bayesian  approach.  Here, 
the  random  parameter  is  described  by  a  known  prior  probability  density  function  (PDF),  and 
the  goal  is  to  estimate  a  particular  realization.  Common  optimization  criteria  for  Bayesian 
estimators  include  minimizing  mean  square  error  (MSE)  between  parameter  and  estimate 
as  well  as  maximizing  a  posteriori  probability.  Both  classical  and  Bayesian  estimators  can 
sometimes  be  difficult  to  implement,  requiring  multidimensional  integration  or  intensive 
iterative  optimization.  In  many  cases,  constraining  the  estimator  to  be  linear  allows  for 
substantial  simplification  and  ease  of  analysis.  Two  well-known  linear  methods  are  least 
squares  and  linear  minimum  MSE  estimation,  also  known  as  the  Wiener  filter  [43]. 

Least  squares  is  a  classical  estimation  method  first  used  by  Gauss  to  study  planetary 
motions  in  1795  [43].  The  goal  is  to  find  the  object  estimate  that  minimizes  the  squared  error 
between  the  given  distorted  image  and  some  deterministic  image  model.  No  probabilistic 
assumptions  are  made  about  the  data  [43].  Actual  performance  is  dependent  on  two  factors: 
the  noise  properties  of  the  distorted  image  and  the  image  model  accuracy.  For  example,  low 
light  images  degraded  by  atmospheric  turbulence  are  not  good  candidates  for  least  squares 
processing  when  the  image  model  or  blur  cannot  be  estimated  accurately.  In  this  case,  a 
wavefront  sensor  (WFS)  can  be  used  to  provide  an  accurate  estimate  of  the  pupil  plane 
phase  aberration.  The  phase  measurement  can  then  be  used  to  estimate  the  blur  function. 
This  technique  is  known  as  deconvolution  from  wavefront  sensing  (DWFS)  [10, 17, 66]  and 
is  based  on  the  least  squares  paradigm.  A  statistical  noise  model  is  not  incorporated  in 
the  traditional  DWFS  estimator  to  suppress  noise  effects.  The  standard  solution  is  to  use 
a  regularization  constant  in  the  estimator  denominator  [20] .  The  regularization  constant  is 
adjusted  by  the  user  based  on  perceived  image  quality. 

The  Wiener  filter  minimizes  MSE  between  the  true  object  and  object  estimate  [43].  It 
was  first  derived  for  two  dimensional  images  by  Helstrom  [29]  and  Slepian  [80].  Their  research 
followed  the  paradigm  established  in  the  seminal  work  by  Norbert  Wiener  with  stationary 
time  series  [94].  Here,  the  object  and  noise  are  assumed  to  be  wide  sense  stationary  random 
processes.  In  this  dissertation,  the  term  stationary  will  refer  to  a  wide  sense  stationary 
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random  process.  A  two  dimensional  stationary  random  process  has  a  constant  mean  and 
autocorrelation  that  is  only  dependent  on  the  distance  between  pixels,  not  the  individual 
pixel  locations  [62].  In  the  Fourier  domain,  a  stationary  random  process  has  uncorrelated 
spatial  frequency  components  [62].  Since  the  Fourier  components  are  uncorrelated,  the 
Wiener  filter  reduces  to  a  simplified  scalar  form  which  weights  each  distorted  image  spatial 
frequency  independently  to  produce  the  estimated  object  spectrum.  However,  many  practical 
image  ensembles  are  non-stationary  [39]  and  have  correlated  Fourier  components  [62].  For 
example,  consider  the  image  of  a  satellite  in  low-earth  orbit  as  collected  by  a  ground-based 
telescope.  Multiple  images  of  the  satellite  are  collected  from  different  perspectives  as  it  passes 
over  the  observation  site.  Thus,  the  mean  object  is  a  blurred  version  of  the  true  satellite 
against  the  black  background  of  space.  Clearly,  this  object  random  process  is  not  constant 
mean.  In  addition,  it  is  not  uncommon  to  use  a  support  constraint  when  processing  these 
images  which  generates  a  non-stationary  image  domain  covariance  [14].  The  scalar  Wiener 
filter  is  not  capable  of  incorporating  complete  object  and  noise  Fourier  domain  correlations. 
Thus,  it  is  sub-optimal  with  respect  to  MSE  for  non-stationary  image  ensembles. 

As  noted  above,  nonlinear  classical  and  Bayesian  estimation  techniques  often  involve 
multidimensional  integration  and  intensive  iterative  optimization.  Linear  techniques  can 
offer  advantages  related  to  computational  savings  and  analysis,  often  at  the  expense  of 
performance.  Thus,  the  desire  to  enhance  the  performance  of  linear  reconstruction  techniques 
is  the  key  motivation  for  solving  the  problem  addressed  in  this  dissertation: 

Develop  enhanced  linear  reconstruction  filters  for  non-stationary  image  ensembles 

incorporating  a  priori  blur  and  noise  models.  Investigate  performance  limitations 

associated  with  imaging  through  atmospheric  turbulence. 

Both  a  classical  and  a  Bayesian  approach  are  addressed  in  this  dissertation.  The  classical 
approach  involves  the  application  of  constrained  least  squares  (CLS)  to  DWFS.  CLS  incor¬ 
porates  a  priori  knowledge  of  the  imaging  scenario  to  constrain  the  set  of  possible  object 
estimates.  The  new  algorithm  incorporates  blur  and  noise  models  to  select  the  appropriate 
regularization  constant  automatically  [11].  No  ad  hoc  regularization  adjustment  is  required. 
CLS  processing  of  DWFS  data  is  demonstrated  using  simulated  satellite  objects  degraded  by 
atmospheric  turbulence  and  measurement  noise.  Measurement  noise  is  defined  as  the  com- 
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bined  contribution  of  shot  noise  and  signal-independent  detector  noise  [14].  The  Bayesian 
approach  involves  the  development  of  a  new  vector  Wiener  filter  which  is  optimal  with  re¬ 
spect  to  MSE  for  a  non-stationary  object  class  degraded  by  atmospheric  turbulence  and 
measurement  noise.  Here,  the  term  vector  alludes  to  the  dependence  of  each  estimated 
Fourier  component  on  other  distorted  image  Fourier  components.  It  should  be  noted  that 
the  general  Wiener  filter  is  the  steady-state  constant-gain  version  of  the  Kalman  filter  for 
shift-invariant  image  models  and  stationary  noise  models  [43].  Thus,  it  is  valid  to  refer  to  the 
new  vector  Wiener  filter  as  a  Kalman  filter  when  the  noise  is  non-stationary.  However,  pre¬ 
vious  image  processing  research  has  included  reference  to  a  Wiener  filter  in  this  application. 
Pratt  proposed  a  generalized  vector  Wiener  filter  that  is  optimal  with  respect  to  MSE  for 
non-stationary  object  ensembles  in  signal-independent  noise  [63].  This  theory  was  extended 
to  images  degraded  by  both  known  blur  and  signal-independent  noise  [65,76].  The  research 
presented  in  this  dissertation  will  use  the  term  Wiener  filter  when  referring  to  this  image 
reconstruction  application.  This  work  involves  the  first  extension  of  this  theory  to  account 
properly  for  shot  noise  [14].  Vector  Wiener  filter  performance  is  also  investigated  when  blur 
statistics  are  substituted  for  exact  knowledge  of  the  blur  function  [12,13]. 

This  dissertation  is  organized  into  seven  chapters.  Chapter  I  presents  the  justification 
for  pursuing  the  suggested  study,  details  the  problem  to  be  solved,  and  outlines  significant 
results.  Chapter  H  contains  background  material  associated  with  atmospheric  turbulence 
and  image  reconstruction.  Chapter  III  outlines  a  new  CLS  processing  technique  for  DWFS. 
Chapter  IV  contains  the  complete  derivation  of  a  new  vector  Wiener  filter  which  incorporates 
model-based  statistical  knowledge  of  object,  blur,  and  noise.  Chapter  V  illustrates  vector 
Wiener  filter  performance  on  astronomical  images  degraded  by  atmospheric  turbulence  and 
noise.  Chapter  VI  presents  vector  Wiener  filter  performance  and  robustness  data  for  gen¬ 
eralized  object  and  blur  models.  Conclusions  and  recommendations  for  further  research  are 
found  in  Chapter  VH.  Mathematical  details  not  included  in  the  main  text  are  compiled  in 
the  Appendices. 
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1.2  Justification  for  Conducting  the  Proposed  Research 

Today,  the  United  States  faces  an  ever-growing  number  of  potential  adversaries  with 
satellite  launch  capability.  Clear,  resolvable  images  of  space  objects  from  ground-based 
telescopes  are  an  absolute  requirement  to  determine  an  opponent’s  intentions  in  space.  Thus, 
the  Air  Force  has  a  requirement  for  high  resolution  imagery  of  earth-orbiting  objects  as  part 
of  its  space  surveillance  mission.  In  general,  two  broad  classes  of  techniques  are  used  to 
increase  Fourier  domain  signal-to-noise  (SNR)  ratio  in  astronomical  images:  (1)  pre-detection 
processing  via  adaptive  optics  (AO)  [27]  and  (2)  post-detection  processing  such  as  speckle 
imaging  [44,47,50]. 

AO  compensates  for  atmospheric  turbulence-induced  wavefront  aberrations  in  real  time 
before  the  light  is  detected  at  the  image  plane.  The  important  components  of  an  AO  system 
are  the  deformable  mirror  (DM),  wavefront  sensor  (WFS),  and  actuator  control  computer 
[74].  Voltages  applied  to  the  DM  actuators  allow  its  figure  to  be  changed  in  real  time. 
The  WFS  senses  the  aberrations  in  the  incoming  wave  by  measuring  gradients  in  small 
subapertures  of  the  telescope  pupil  [91].  This  information  is  then  sent  to  the  actuator  control 
computer  which  adjusts  the  DM  to  apply  an  estimate  of  the  conjugate  of  the  wavefront 
aberration.  The  correction  imposed  by  the  DM  cancels  out  the  aberration,  leading  to  a 
narrower  blur  or  point  spread  function  (PSF)  and  an  improved  image.  This  process  must 
occur  at  speeds  on  the  order  of  the  rate  of  change  of  the  wavefront  aberration  to  be  effective 
[91].  Typically,  these  speeds  range  from  approximately  tens  of  Hertz  to  a  few  hundred  Hertz 
[91].  Figure  1.1  illustrates  the  installation  of  an  AO  system  on  a  ground-based  telescope. 

The  first  work  to  address  the  post-detection  processing  of  images  degraded  by  atmo¬ 
spheric  turbulence  was  speckle  imaging  [44,47,50].  The  term  speckle  refers  to  the  data, 
which  consist  of  a  set  of  short  exposure,  speckled  images.  In  this  context,  short  exposure 
refers  to  a  sufficiently  short  integration  time  to  freeze  an  individual  realization  of  the  at¬ 
mospheric  turbulence-induced  aberration  in  the  image  measurement.  In  this  technique,  the 
object  is  usually  estimated  by  first  estimating  the  modulus  and  phase  of  its  Fourier  trans¬ 
form.  Labeyrie  showed  that  the  squared  modulus  of  the  object  Fourier  transform  could  be 
estimated  from  a  large  set  of  short  exposure  images  [47].  The  method  requires  an  ensemble 
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Figure  1.1 


Diagram  of  a  typical  AO  system  as  part  of  a  large  telescope. 
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Figure  1.2  Block  diagram  of  the  speckle  imaging  process. 

of  reference  point  source  images  along  with  the  data  set.  The  reference  images  are  used  to 
estimate  the  Fourier  transform  of  the  atmospheric-optical  system  PSF  or  optical  transfer 
function  (OTF).  The  squared  modulus  of  both  data  and  OTF  estimates  are  averaged  to 
reduce  noise.  Then,  the  average  squared  modulus  of  the  OTF  estimates  is  used  in  a  decon¬ 
volution  procedure  to  estimate  the  squared  modulus  of  the  object.  This  estimation  scheme  is 
possible  because  the  squared  modulus  of  the  OTF  is  non-zero  out  to  the  diffraction-limited 
cutoff  of  the  optical  system  [74].  While  Fourier  modulus  information  can  be  extracted  us¬ 
ing  the  above  technique,  phase  information  is  usually  required  to  form  a  usable  image  [74]. 
Two  methods  are  commonly  used  to  extract  the  Fourier  phase  from  the  data  ensemble:  the 
Knox-Thompson  [44]  and  bispectrum  [50]  methods.  Both  methods  are  based  on  the  fact  that 
certain  higher  order  moments  of  the  complex  Fourier  transform  of  speckled  images  contain 
encoded  information  about  the  object  phase.  Figure  1.2  gives  a  block  diagram  of  the  speckle 
imaging  process. 
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Figure  1.3  Block  diagram  of  DWFS. 


A  third  class  of  imaging  techniques  involves  the  combination  of  both  AO  and  post¬ 
detection  processing  and  is  known  as  hybrid  imaging  [74].  A  notable  hybrid  technique  is 
DWFS  [10,17,66].  In  DWFS,  a  WFS  is  used  to  measure  the  pupil  plane  phase  aberration 
associated  with  each  short  exposure  image.  The  phase  measurement  is  used  to  estimate  the 
OTF.  The  short  exposure  image  and  OTF  estimate  can  then  be  used  to  estimate  the  ob¬ 
ject  via  a  deconvolution  filter.  The  DWFS  technique  was  first  proposed  by  Fontanella  [10], 
extended  by  Pried  [17],  and  further  developed  by  Primot  et  al.  [66].  Primot  et  al.  also 
conducted  the  first  performance  analysis  of  DWFS  [66].  A  variant  of  their  estimator  was 
later  validated  on  astronomical  data  [20].  Welsh  and  Von  Niederhausern  further  investigated 
DWFS  performance  by  incorporating  an  optimal  wavefront  phase  estimate  [93].  Roggemann 
et  al.  [75]  showed  that  the  Primot  estimator  was  biased  and  suggested  a  related  unbiased 
estimator.  Roggemann  and  Welsh  also  derived  an  SNR  expression  for  DWFS  [73]  and  con¬ 
ducted  further  comparison  with  speckle  imaging  and  traditional  linear  deconvolution  [92]. 
Figure  1.3  gives  a  block  diagram  of  the  DWFS  process. 
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Speckle  imaging  techniques  generally  incorporate  a  deconvolution  filter  to  estimate  the 
object  modulus  [74].  AO  compensated  images  also  benefit  from  linear  reconstruction  to 
deconvolve  blurring  due  to  the  attenuation  of  high  spatial  frequencies  in  the  compensated 
images  [70,72].  Typical  linear  reconstruction  methods  such  as  the  inverse  and  pseudo- Wiener 
filters  [21]  require  an  explicit  estimate  of  the  OTF.  A  priori  knowledge  of  the  object  class  and 
noise  is  typically  not  used  to  deconvolve  compensated  images.  CLS  processing  of  DWFS  data 
and  the  vector  Wiener  filter  offer  the  potential  to  use  statistical  model-based  information 
about  the  OTF  and  noise  to  improve  images  from  Air  Force  ground-based  surveillance  sites. 
As  of  this  writing,  no  sources  have  been  found  in  the  literature  which  document  the  use  of 
CLS  estimation  to  process  DWFS  data.  However,  Primot  et  al.  note  the  potential  of  such 
a  scheme  in  their  key  paper  [66].  Similarly,  no  sources  are  available  which  document  the 
use  of  a  Fourier  domain  linear  minimum  MSE  estimator  to  process  non-stationary  image 
ensembles  degraded  by  random  blur  and  shot  noise.  Several  researchers  have  addressed 
the  problem  of  image  reconstruction  in  the  presence  of  random  blur  but  only  in  the  image 
domain  [4,25,26,87].  The  potential  benefits  of  statistical  model-based  PSF  information  has 
been  mentioned  by  prominent  researchers  in  the  area  of  blind  deconvolution.  Schulz  [77] 
has  recommended  the  use  of  PSF  statistical  models  in  his  algorithms  when  this  information 
is  available.  However,  his  work  and  the  work  of  other  researchers  in  this  field  have  not 
considered  the  degrading  PSF  as  a  random  quantity  [31,33,48,77,81].  Thus,  the  research 
outlined  in  this  dissertation  makes  a  unique  contribution  to  a  critical  Air  Force  mission  and 
the  field  of  image  reconstruction. 

1.3  Approach 

The  problem  statement  is  addressed  in  two  ways.  This  two-prong  approach  is  consistent 
with  the  natural  classical-Bayesian  division  prevalent  in  statistical  estimation  theory  [43]. 
First,  the  application  of  CLS  to  DWFS  is  investigated.  CLS  represents  a  classical  estimation 
approach  since  the  object  is  assumed  to  be  a  deterministic,  yet  unknown  quantity.  Only  the 
constraint  incorporates  model-based  statistical  information.  Second,  a  new  linear  Bayesian 
estimator,  referred  to  as  the  vector  Wiener  filter,  is  derived.  In  this  case,  model-based 
statistical  knowledge  of  object,  blur,  and  noise  is  assumed. 
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The  traditional  DWFS  deconvolution  filter  given  by  Primot  et  al.  [66]  provides  un¬ 
acceptable  noise  amplification  when  processing  low  light  image  ensembles.  The  standard 
solution  is  to  add  a  parametric  regularization  constant  [20]  or  SNR  term  [75]  to  the  filter 
denominator.  These  approaches  are  analogous  to  CLS  [37]  and  the  parametric  Wiener  fil¬ 
ter  [21],  respectively.  The  regularization  is  adjusted  by  the  user  based  on  the  perceived  image 
quality.  Thus,  the  resultant  object  estimate  does  not  satisfy  any  mathematical  optimality 
criterion.  In  this  dissertation,  a  modified  CLS  estimation  scheme  is  developed  which  provides 
optimal  processing  of  noisy  DWFS  data  [11].  Here,  optimal  refers  to  the  object  estimate 
which  minimizes  a  CLS  objective  function  incorporating  DWFS  data.  Unlike  previous  CLS 
algorithms  [37, 71] ,  this  approach  incorporates  ensemble  average  data  directly  to  reduce  noise 
effects.  The  solution  uses  the  Lagrange  multiplier  technique  [21,37].  A  closed  form  solution 
for  the  object  estimate  is  obtained  which  is  analogous  to  the  traditional  DWFS  deconvolution 
filter  [20,75]  with  the  Lagrange  multiplier  serving  as  a  regularization  constant.  An  iterative 
approach  based  on  Newton-Raphson  minimization  [40,  71]  is  used  to  find  the  appropriate 
regularization  constant.  The  iteration  incorporates  the  statistics  of  both  the  OTF  and  noise. 

CLS  processing  of  noisy  DWFS  data  relies  on  WFS  hardware  and  iterative  processing 
to  deconvolve  turbulence  effects.  A  non-iterative  Bayesian  approach  to  the  reconstruction 
of  astronomical  images  is  also  developed,  known  as  the  vector  Wiener  filter.  As  noted 
previously,  a  scalar  filter  weights  each  Fourier  component  of  the  distorted  data  independently. 
In  contrast,  a  vector  filter  incorporates  many  Fourier  components  of  the  distorted  image 
to  estimate  a  given  Fourier  component  of  the  object.  The  appropriate  weighting  of  each 
component  is  determined  by  the  object,  OTF,  and  noise  correlations.  First,  a  vector  Wiener 
filter  is  derived  that  properly  accounts  for  a  random  OTF  and  shot  noise  effects  [14].  This 
new  linear  filter  has  the  advantage  of  incorporating  object,  OTF,  and  shot  noise  correlations 
between  different  spatial  frequency  components.  When  a  scalar  Wiener  filter  is  applied  to 
a  non-stationary  image  ensemble,  this  correlation  information  is  not  used.  The  result  is  a 
sub-optimal  solution  with  respect  to  MSE  when  compared  to  the  vector  Wiener  filter.  Next, 
vector  Wiener  filter  performance  is  investigated  for  both  a  fixed  and  random  OTF  [12-14]. 
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Finally,  filter  performance  and  robustness  are  examined  for  generalized  object  and  blur 
models. 

This  research  effort  relies  on  theoretical  analysis  and  Monte  Carlo  simulation.  The 
Monte  Carlo  simulation  for  the  atmospheric  turbulence-degraded  images  is  based  on  a 
Fourier-series  based  phase  screen  generator  developed  by  Dr.  Byron  Welsh  [89].  This  new 
phase  screen  generator  properly  models  the  spatial  and  temporal  correlations  between  wave- 
front  phase  screens  based  on  von  Karman  statistics  [89].  However,  this  study  is  only  con¬ 
cerned  with  the  spatial  correlation  of  the  OTF.  Thus,  the  temporal  capability  of  the  phase 
screen  generator  is  not  used.  Performance  comparison  between  the  scalar  and  vector  Wiener 
filters  is  an  important  part  of  this  study.  Here,  estimator  performance  is  based  on  visual 
image  comparison,  MSE,  and  mean  square  phase  error  (MSPE).  MSE  at  a  given  image  pixel 
is  defined  as  the  expected  value  of  the  squared  difference  between  the  true  object  and  the 
object  estimate  at  that  pixel.  Similarly,  MSPE  for  a  given  Fourier  component  is  the  expected 
value  of  the  squared  difference  between  the  true  object  phase  and  the  estimated  object  phase 
for  that  Fourier  component. 

1.4  Scope  and  Assumptions 

Two  new  applications  of  model-based  statistical  knowledge  to  linear  filter  theory  are 
presented  in  this  dissertation.  The  primary  study  variables  are  the  object  class,  light  level, 
detector  read  noise  variance,  and  turbulence  strength.  The  emphasis  here  is  on  analysis  and 
simulation. 

In  the  discussion  presented  in  Section  1.3,  the  object  irradiance  distribution  was  as¬ 
sumed  to  be  a  random  process  with  known  spatial  frequency  statistics.  The  concept  of  a 
random  object  in  image  reconstruction  is  not  new  [29, 63, 80]  and  is  critical  to  a  Bayesian 
development.  In  Chapters  IV  and  V,  perfect  a  priori  knowledge  of  the  object  class  statistics 
is  assumed.  For  example,  one  common  class  of  astronomical  objects  is  the  binary  star  pair. 
A  priori  knowledge  could  include  the  number  of  components  (two),  ratio  between  primary 
and  secondary  component  irradiance,  and  object  support.  In  this  case,  exact  knowledge  of 
the  true  object  irradiance  distribution  is  unavailable  since  the  filter  has  no  knowledge  of  the 
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component  separation  or  orientation.  In  general,  the  statistical  object  model  can  be  viewed 
as  a  constraint  on  the  filter  output  [43].  Clearly,  the  object  statistics  can  never  be  known 
perfectly  in  a  real  application.  Thus,  filter  robustness  is  studied  in  Chapter  VI. 

1.5  Significant  Contributions  and  Results 

This  section  highlights  overall  contributions  and  results  associated  with  this  disserta¬ 
tion  research.  Here,  the  logical  questions  “What  is  new?”  and  “What  is  important?”  are 
answered.  The  respective  chapters  are  listed  to  aid  the  reader  in  finding  topics  of  interest. 

•  Chapter  III  outlines  a  new  application  of  CLS  to  the  DWFS  processing  of  low  light 
images.  The  technique  is  practical  and  computationally  inexpensive.  In  all  cases  exam¬ 
ined,  the  Newton-Raphson  iteration  converged  to  a  solution  in  less  than  10  iterations. 

•  Chapter  IV  presents  the  derivation  of  a  new  vector  Wiener  filter  incorporating  the 
semi-classical  model  of  photoelectric  light  detection.  The  filter  uses  complete  OTF 
and  shot  noise  statistical  models. 

•  Chapter  V  presents  the  first  application  of  second  order  OTF  statistics  between  dif¬ 
ferent  spatial  frequencies  in  a  Wiener  filter.  These  OTF  statistics  are  associated  with 
imaging  through  atmospheric  turbulence. 

•  Chapter  VI  contains  the  first  performance  study  to  establish  quantitative  limits  on  the 
vector  Wiener  filter  associated  with  object  and  OTF  statistical  models.  The  object 
spatial  SNR  and  correlation  coefficient  provide  a  fundamental  limit  on  vector  Wiener 
filter  MSE  performance.  For  some  object  classes,  the  OTF  SNR  also  limits  MSE 
performance. 

•  Chapter  VI  contains  the  first  robustness  study  of  the  vector  Wiener  filter  with  respect 
to  object  model  error.  Error  in  the  object  spatial  SNR  produces  a  substantial  increase 
in  MSE. 
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1.6  Summary 

In  this  chapter,  the  dissertation  problem  has  been  presented,  which  is  to  develop  linear 
filters  for  non-stationary  image  ensembles  incorporating  blur  and  noise  statistical  models. 
This  research  is  justified  based  on  the  Air  Force  requirement  for  clear,  resolvable  images  of 
space-based  objects.  Linear  reconstruction  provides  an  important  function  by  deconvolving 
images  processed  via  more  sophisticated  techniques  and  hardware.  The  approach  outlined 
in  this  dissertation  is  to  investigate  both  CLS  processing  of  DWFS  data  and  a  new  vector 
Wiener  filter.  Both  techniques  incorporate  blur  and  noise  statistics.  In  addition,  the  vector 
Wiener  filter  relies  on  model-based  statistical  knowledge  about  the  object  class.  Before 
presenting  these  techniques  in  Chapters  III  and  IV,  the  next  chapter  provides  important 
background  regarding  atmospheric  turbulence  and  key  image  reconstruction  methods. 


1-13 


IL  Background 


2.1  Introduction 

The  primary  objective  of  this  chapter  is  to  review  background  material  related  to  gen¬ 
eral  image  reconstruction  concepts  and,  more  specifically,  the  reconstruction  of  astronomical 
images.  The  image  degradation  model  associated  with  a  noisy,  turbulence-degraded  image 
is  introduced  as  well  as  information  about  atmospheric  turbulence.  Unless  otherwise  noted, 
this  image  model  will  be  used  throughout  the  dissertation.  The  rest  of  the  chapter  will 
briefly  review  some  important  classical  and  Bayesian  techniques  in  image  reconstruction.  In 
each  case,  the  strengths  and  weaknesses  of  past  research  will  be  highlighted. 

The  rest  of  this  chapter  is  organized  as  follows.  Section  2.2  presents  the  image  degrada¬ 
tion  model.  Atmospheric  turbulence  theory  is  reviewed  in  Section  2.3  with  emphasis  on  ex¬ 
pressions  for  the  optical  transfer  function  (OTF)  statistics.  Classical  estimation  schemes  are 
addressed  in  Section  2.4  to  include  least  squares,  maximum  likelihood  (ML),  and  Gerchberg- 
Saxton  iteration.  Section  2.5  introduces  Bayesian  estimation  concepts  to  include  the  Wiener 
filter,  Kalman  filter,  and  maximum  a  posteriori  (MAP)  estimation.  Section  2.6  provides 
a  brief  summary  of  the  chapter  and  relates  this  past  research  to  the  work  outlined  in  this 
dissertation. 

2.2  Image  Degradation  Model 

The  standard  model  for  a  noise-free  linear  shift-invariant  imaging  system  can  be  written 
as  [22] 

i{x, y)  =  h{x, y)  *  o(x,  y),  (2.1) 

where  i{x,y)  is  the  noiseless  image,  h{x,y)  is  the  impulse  response  or  point  spread  function 
(PSF)  associated  with  the  blur,  o(x,y)  is  the  object,  {x,y)  is  a  discrete  point  in  the  image 
domain,  and  *  denotes  convolution.  Taking  the  Fourier  transform  of  both  sides  of  Eq.  (2.1) 
yields 

I{u,v)  =  7i{u,v)0{u,v),  (2.2) 
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where  a  capital  letter  denotes  the  discrete  spatial  Fourier  transform  of  an  associated  lower 
case  quantity  and  W(u,  v)  is  the  Fourier  transform  of  h(x,  y)  also  known  as  the  OTF.  The 
point  (u,  v)  denotes  a  discrete  spatial  frequency. 

The  OTF  for  a  diffraction-limited  full  aperture  is  a  deterministic,  tapered  low-pass  fil¬ 
ter.  Random  aberrations  in  the  optical  system  pupil  due  to  atmospheric  turbulence  further 
degrade  performance,  resulting  in  a  more  attenuated  OTF,  especially  at  high  spatial  frequen¬ 
cies  [22].  To  complicate  matters,  the  noiseless  image  i{x,y)  is  typically  not  available.  An 
image  model  must  account  for  additional  degradation  due  to  the  detection  process.  Photon- 
matter  interactions  in  light  detectors  are  random  and  require  a  statistical  description.  The 
semi-classical  model  is  based  on  three  assumptions  about  photon  statistics  [74]: 

1.  The  probability  of  occurrence  of  a  single  photoevent  in  a  small  area  dA  during  a  time 
interval  dt  is 

P{l,dt,dA)  =  T)  dt  dAi(x,y,t),  (2.3) 

where  i{x,y,t)  =  h{x,y,t)  *  o{x^y^t),  dA  is  small  compared  to  the  coherence  area  of 
the  light,  dt  is  short  compared  to  the  coherence  time  of  the  light,  and  77  is  the  quantum 
efficiency  of  the  detector. 

2.  The  probability  of  more  than  one  photoevent  occurring  in  the  area  dA  during  time 
interval  dt  is  very  small  compared  to  the  probability  associated  with  either  one  or  zero 
photoevents. 

3.  The  number  of  photoevents  K  occurring  in  non-overlapping  space  or  time  intervals  is 
statistically  independent. 

Based  on  these  assumptions,  the  random  variable  K  is  governed  by  Poisson  statistics  [23,74]. 

Poisson  random  process  sample  functions  consist  of  sets  of  Dirac  delta  functions  [62] . 
Thus,  a  photon-limited  image  is  defined  as  [74] 

K 

d{x,y)  =  ^5{x-Xn,y-yn),  (2.4) 

n=l 
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where  {xn,yn)  is  the  location  of  the  photoevent  in  the  image  plane  and  K  is  the  total 
number  of  photoevents  making  up  the  image.  Randomness  is  associated  with  the  number 
and  location  of  the  photoevents.  The  randomness  considered  by  the  semi-classical  model 
is  referred  to  as  photon  or  shot  noise  and  is  a  form  of  signal-dependent  noise.  Signal- 
dependent  refers  to  the  situation  in  which  the  strength  of  the  noise  depends  on  the  number 
and  distribution  of  photoevents  [74].  Photon  noise  typically  imposes  more  severe  limitations 
than  diffraction,  especially  at  low  light  levels.  The  Dirac  delta  functions  of  Eq.  (2.4)  are 
discontinuous  in  the  image  domain,  which  presents  difficulties  when  conducting  statistical 
analysis.  Thus,  we  are  motivated  to  analyze  imaging  performance  in  the  Fourier  domain  via 
the  Fourier  transform  of  Eq.  (2.4)  which  can  be  written  as 

K 

D{u,  ?;)  =  XI  {-i27r(«a;„  -f-  vyn)}  .  (2.5) 

n=l 

Signal-independent  noise  is  also  present  in  many  detectors  used  for  image  collection.  For 
example,  charge-coupled  device  (CCD)  detector  output  is  degraded  by  signal-independent 
additive  noise  known  as  read  noise  [74].  In  this  dissertation,  signal-independent  noise  will 
be  represented  by  the  random  variable  Up  having  the  following  statistical  properties: 

1.  Up  is  zero  mean. 

2.  rip  is  spatially  uncorrelated  with  uniform  variance 

3.  Up  is  statistically  independent  of  K  and  {xn,yn)- 

An  image  model  which  properly  accounts  for  both  signal-dependent  and  signal-independent 
noise  can  now  be  written  as  [74] 

K  p 

d{x,y)  =  ^5{x  -  Xn,y  -  yn)  +  ~  ~  Vp)^  (2-6) 

n=l  p=l 

where  {xp,  yp)  is  the  location  of  the  image  pixel,  P  is  the  total  number  of  pixels  in  the 
detector  array,  and  d{x,  y)  now  represents  a  detected  image  as  collected  by  a  CCD  detector. 
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The  corresponding  Fourier  domain  expression  can  be  written  as 


K  p 

D{u,  v)  =  '^  exp  {-j2Ti:{uXn  +  vyn)}  +  Wp  exp  {-j2TT{uXp  +  vyp)}  .  (2.7) 

n=l  p=l 

Equation  (2.7)  is  a  key  expression  used  to  model  noise  effects  in  this  dissertation. 

To  investigate  both  constrained  least  squares  (CLS)  processing  of  deconvolution  from 
wavefront  sensing  (DWFS)  data  and  the  vector  Wiener  filter,  a  vector-matrix  expression 
is  needed  for  Eq.  (2.6).  The  new  expression  includes  an  additive  noise  vector  n  which 
incorporates  both  shot  and  read  noise  effects  based  on  writing  the  total  noise  as  the  difference 
between  the  detected  image  and  the  image  degraded  by  the  PSF  only,  denoted  n{x,  y)  = 
d{x,y)  —  i(x^y)  [14,37].  Thus,  Eq.  (2.6)  can  now  be  written  as  [14] 

d  =:  Ho  +  n,  (2.8) 

where  d  and  o  are  P-length  vector  versions  of  the  functions  d(x,  y)  and  o{x,  y),  respectively. 
The  matrix  if  is  a  P  x  P  block-circulant  matrix  representing  the  shift-invariant  PSF.  H  and 
o  are  properly  ordered  to  perform  the  discrete  convolution  of  Eq.  (2.1)  [21,37].  Equation 
(2.7)  can  also  be  written  using  this  vector-matrix  notation  such  that 

D  =  'HQO  +  TSl,  (2.9) 

where  D,  O,  and  N  are  P-length  vector  versions  of  the  Fourier  domain  functions  D{u,v), 
0{u,v),  and  N{u,v),  respectively.  The  notation  ©  represents  an  entrywise  or  Hadamard 
product  [35]  and  H  is  a  P-length  vector  containing  the  OTF.  Equations  (2.8)  and  (2.9)  are 
key  expressions  in  the  development  of  CLS  processing  of  DWFS  data  and  the  vector  Wiener 
filter  in  Chapters  III  and  IV,  respectively. 

The  next  section  reviews  atmospheric  turbulence  theory.  Expressions  for  the  mean  and 
spatial  correlation  of  the  turbulence-degraded  OTF  are  presented.  These  OTF  statistics  will 
play  an  important  role  in  both  CLS  processing  of  DWFS  data  and  the  vector  Wiener  filter. 
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2.3  Atmospheric  Turbulence 


An  undergraduate  Physics  textbook  provides  the  necessary  tools  to  predict  the  resolu¬ 
tion  of  an  imaging  system.  Theoretically,  the  minimum  resolvable  angle  seen  by  a  telescope 
is  limited  by  diffraction  and  can  be  expressed  as 


_  1.22A 

“  D  ’ 


(2.10) 


where  9  is  the  minimum  resolvable  angle  in  degrees,  A  is  the  mean  wavelength  of  the  incident 
light,  and  D  is  the  diameter  of  the  aperture.  The  over  line  notation  denotes  the  expectation 
operator  applied  to  the  designated  quantity.  As  noted  in  the  previous  section,  ground-based 
imaging  systems  rarely  achieve  diffraction-limited  performance.  Instead,  the  atmosphere 
imposes  a  fundamental  limit  on  spatial  resolution.  Atmospheric  turbulence  affects  imaging 
systems  by  causing  both  spatial  and  temporal  random  fluctuations  in  the  index  of  refraction 
of  the  atmosphere.  These  index  of  refraction  fluctuations  impose  random  phase  aberrations 
on  the  incoming  light  [23,69].  The  primary  consequence  of  these  random  phase  aberrations 
is  a  general  broadening  of  the  PSF  which  manifests  itself  as  blurring  and  lowered  resolution 
when  compared  to  the  system  predicted  by  Eq.  (2.10). 

Atmospheric  turbulence  is  caused  by  turbulent  air  motion.  The  source  of  this  air 
motion  is  the  heating  and  cooling  of  the  Earth  by  the  sun.  Large  air  masses  gain  heat 
directly  from  the  sun  during  the  day.  At  night,  heat  is  also  coupled  to  these  air  masses 
as  the  Earth  cools.  As  a  result,  large-scale  temperature  variations  are  produced.  These 
temperature  variations  lead  to  pressure  differences  which  result  in  large  scale  air  motion. 
Initial  large  scale  air  motions  break  down  into  smaller  and  smaller  scale  motions  until  the 
atmosphere  is  distributed  into  randomly  sized  pockets  of  air,  each  with  its  own  temperature, 
as  shown  in  Fig.  2.1.  These  pockets  of  air  are  called  turbulent  eddies  [69].  Since  the  index 
of  refraction  of  air  is  dependent  on  temperature,  the  atmosphere  has  a  non-uniform  index  of 
refraction. 


2.3.1  Turbulence  Statistics.  Atmospheric  turbulence  creates  a  medium  which  has  a 
non-uniform  or  random  index  of  refraction  associated  with  the  distribution  of  the  turbulent 
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Figure  2.1  Turbulent  eddies  and  their  effects  on  an  incoming  unperturbed  plane  wave. 

eddies.  Thus,  statistical  models  are  required  to  understand  turbulence  effects  fully.  Tur¬ 
bulence  modeling  has  received  much  attention  in  the  literature.  However,  most  published 
research  flows  from  the  key  results  by  Kolmogorov  [45],  Tatarskii  [83],  and  Pried  [15,16].  Kol¬ 
mogorov  provided  a  statistical  model  related  to  spatial  structure  in  turbulent  air  flows  [45]. 
Tatarskii  used  Kolmogorov’s  results  to  model  wave  propagations  through  random  index  of 
refraction  distributions  [83].  Pried  applied  and  extended  Tatarskii’s  work  to  optical  propa¬ 
gation  problems  [15, 16]. 

Kolmogorov  theory  gives  a  mathematical  description  for  index  of  refraction  fluctuations 
[45].  The  index  of  refraction  spatial  power  spectral  density  (PSD)  provides  a  frequency 
domain  statistical  model  for  the  number  and  size  of  the  turbulent  eddies  and  can  be  written 
as  [74] 

«^M=0.033CJ(^)k-^,  (2.11) 
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where  the  superscript  K  denotes  the  Kolmogorov  spectrum,  k  is  the  scalar  wavenumber, 
is  a  measure  of  turbulence  strength  called  the  structure  constant,  and  the  subscript  n  denotes 
the  atmospheric  index  of  refraction  [23] .  The  constants  Lg  and  Ig  denote  the  outer  and  inner 
scale,  respectively.  These  quantities  represent  the  characteristic  dimension  of  the  largest 
and  smallest  turbulent  eddies  [74].  In  many  practical  systems,  the  turbulence  strength  is  a 
function  of  the  distance  from  the  imaging  aperture,  denoted  by  z.  Equation  (2.11)  is  not 
a  useful  model  for  the  index  of  refraction  PSD  when  0  because  of  the  non-integrable 
singularity  at  k  =  0.  An  alternate  form  known  as  the  von  Karman  spectrum  is  finite  for  all 
K>0  and  can  be  expressed  as  [74] 


0.033C;(z) 

+ ^2)11/6  > 


(2.12) 


where  Kg  =  2'K/Lg  and  the  superscript  V  represents  the  von  Karman  spectrum.  The  von 
Karman  statistical  model  will  be  used  in  theoretical  development  and  simulation  throughout 
this  dissertation. 

Equation  (2.12)  is  important  in  deriving  statistical  models  for  the  turbulence-induced 
perturbations  on  a  propagating  wavefront.  Two  limiting  cases  are  commonly  studied:  near 
field  and  far  field.  In  the  near  field  case,  only  the  perturbations  affecting  the  wavefront 
phase  are  considered.  Under  the  far  field  assumption,  both  amplitude  and  phase  effects 
are  modeled.  Here,  it  is  assumed  that  the  wavefront  amplitude  perturbations  are  negligible 
compared  to  the  phase  perturbations.  This  assumption  is  commonly  used  in  a  standard 
geometrical  optics  model  [74].  Thus,  only  the  near  field  case  is  considered  in  this  dissertation. 
Using  a  layered  turbulence  model  and  assuming  the  index  of  refraction  fluctuations  are  a 
Gaussian  random  process  [15,23,74],  it  is  possible  to  derive  an  expression  for  the  spatial 
correlation  function  of  the  phase  perturbations.  Consider  an  incident  wavefront  and  let 
(j){x,y)  indicate  the  wavefront  phase  perturbations  in  the  optical  system  pupil.  Then  the 
spatial  correlation  of  (f>{x,y),  denoted  B^^{Ax,Ay),  is  [74] 
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R^^(x,y]x',y') 


Ay) 


=  ^[(l>{x,y)<l>ix',y')] 


3.089(27r)-®/® 

(L,^{Ax)^  +  {Ay)^ 

25/6r[(ll/6)] 

V 

X-K’s/e 

(2‘KsJ{Axy  +  {AyY\ 

1  ) 

5 

(2.13) 


where  E[«]  denotes  the  expectation  operator,  Aa;  =  x  —  x',  Ay  =  y  —  y',  {x,y)  and  {x',y') 
are  discrete  points  in  the  pupil  plane,  r[(«)]  is  the  Gamma  function,  K^/^[{*)]  is  the  Bessel 
function  of  the  second  kind  of  order  5/6,  and  Tq  is  the  Fried’s  parameter  defined  as  [74] 


To  =  0.185 


(2.14) 


The  Pried  parameter  can  be  interpreted  as  the  seeing  cell  or  aperture  size  beyond  which 
further  increases  in  optical  system  diameter  result  in  no  further  increase  in  resolution  [74]. 
The  spatial  correlation  function  given  in  Eq.  (2.13)  is  the  key  statistical  quantity  used  to 
model  the  effect  of  atmospheric  turbulence  on  imaging  system  performance. 

While  knowledge  of  the  phase  correlation  function  greatly  enhances  our  understanding 
of  turbulence  effects,  a  related  quantity  known  as  the  phase  structure  function  is  of  interest 
in  many  applications.  The  phase  structure  function,  denoted  D^^{Ax,  Ay),  is  defined  as  [23] 


D^^{Ax,  Ay)  =  2R^^{0, 0)  -  2R^^{Ax,  Ay).  (2.15) 


Substituting  the  von  Karman  phase  correlation  function  of  Eq.  (2.13)  into  Eq.  (2.15)  yields 
[74] 
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Ay) 


As  shown  below,  Eq.  (2.16)  is  an  important  quantity  in  the  derivation  of  first  and  second 
order  OTF  statistics. 


2.3.2  Optical  Transfer  Function  Statistics.  In  Section  2.2,  the  semi-classical  model 
for  light  detection  was  introduced,  as  shown  in  Eq.  (2.4).  This  model  incorporates  Dirac 
delta  functions  which  are  discontinuous  in  the  image  domain.  Therefore,  statistical  analysis 
is  more  straightforward  in  the  Fourier  domain.  To  incorporate  model-based  information 
about  the  blur  or  optical  system  PSF  in  this  alternate  domain,  OTF  statistical  expressions 
are  needed.  The  random  OTF  H{u,v)  is  defined  as  [74] 


'H{u,v)  = 


W {uXf,  vXf)  -k  W {uXf,  vXf) 
Nf 


(2.17) 


where  f  is  the  optical  system  focal  length,  u  =  x/(Xf),  v  =  y/{Xf),  Np  =  TF'(0,0):*-TF(0, 0), 
★  denotes  the  two  dimensional  correlation  operation,  and  W{x,y)  is  the  generalized  pupil 
function  which  incorporates  the  phase  aberrations  such  that  [74] 


W{x,y)  =  TFp(x,y)exp{j(^(a;,y)}. 


(2.18) 


Wp{x.iy)  is  a  real-valued  function  describing  the  unaberrated  pupil. 

The  mean  OTF  can  be  derived  using  two  distinct  approaches  [74].  The  first  approach 
relies  on  an  interferometric  view  of  imaging  and  the  Van-Cittert-Zernike  Theorem  [23].  The 
second  uses  a  thin  screen  model  and  proceeds  directly  from  Eq.  (2.17).  Regardless  of  the 
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approach  taken,  the  result  is  the  straightforward  expression  [74] 


Hiu,  v)  =  Ha{u,  v)Ho{u,  v),  (2.19) 

where  T-Ca{u,  v)  is  the  transfer  function  due  to  the  atmospheric  turbulence  and  Ho{u,  v)  is  the 
transfer  function  associated  with  the  aberration-free  optics.  Pried  derived  both  long  exposure 
and  short  exposure  expressions  for  'Ha{u,v)  [15].  The  term  “long  exposure”  refers  to  the 
situation  in  which  the  imaging  system  has  been  exposed  to  many  independent  realizations 
of  the  atmospheric  turbulence.  Here,  the  term  “short  exposure”  refers  to  the  case  in  which 
wavefront  tilt  has  been  removed.  Fried’s  long  exposure  derivation  relies  on  the  assumption 
that  the  phase  perturbation  ^{x,y)  is  a  stationary  Gaussian  random  process  [74].  After  a 
change  of  variables  and  simplification,  the  final  expression  for  the  long  exposure  OTP  is  [74] 

=  exp|-iD^^(uA/,uA/)| .  (2.20) 

The  short  exposure  or  wavefront  tilt-removed  transfer  function  HasEi'^i  is  also  of  interest. 
Here,  the  residual  phase  after  tilt  removal  can  be  written  as  [74] 

(a;,  y)  =  (t>{x,  y)  -  {a^x  +  Oj,y),  (2.21) 

where  ax  and  Oj,  are  coefficients  describing  the  tilt  of  the  wavefront  phase  over  the  pupil.  In  his 
derivation  of  Pried  assumed  that  (f>r{x,y)  was  uncorrelated  with  tilt  coefficients 

ax  and  Oj,  [15].  While  this  assumption  is  not  strictly  valid  from  a  mathematical  viewpoint, 
the  correlation  is  small  when  D/vo  is  large  [74].  Thus,  the  short  exposure  OTP  can  be 
written  as 

^as£,(w,'y)  =  exp  (D^^{uXf,vXf)  -  ^X'^f  (a2  +  (y  +  | .  (2.22) 

Equations  (2.20)  and  (2.22)  provide  theoretical  expressions  for  the  first  order  OTP  statis¬ 
tics.  The  mean  OTP  was  applied  to  the  reconstruction  of  atmospheric  turbulence-degraded 
images  as  early  as  1967  [80].  Second  order  OTP  statistics  have  also  been  applied  to  such 
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problems  as  speckle  imaging  [23,47,74].  In  this  application,  the  speckle  transfer  function 
E{|7Y(m,  is  an  important  theoretical  quantity.  However,  the  value  of  a  priori  knowledge 
of  the  correlations  between  different  OTF  spatial  frequencies  has  yet  to  be  explored.  In  this 
case,  we  require  the  complete  second  order  statistical  quantity 

Rnn{u,  V,  u\  v')  =  E  \H{u,  v)H*{u',  v')] ,  (2.23) 

where  the  superscript  *  denotes  a  complex  conjugate.  Following  the  same  technique  used  to 
derive  an  expression  for  the  speckle  transfer  function  [74],  the  OTF  correlation  function  is 

exp  +  u2)  +  D^^(XfVu'^  + 

Rnniu,v,u,v)  =  - ^ ^ ^ 

^  /  /  /  /  y  ~ 

xWp{x',y')Wpix'  -  u'Xf,y'  -  v'Xf) 

X  exp  {-i  ^£>00  {x  -  x'y  +  {y-  y')^) 

{x-x'  -  uXf  +  u'XfY  -\-{y-y'  -  vXf  +  v'XfY^ 
(x-x'  -  uXfY  +  {y-y'  -  vXfy^ 

{x-x'  +  u'XfY  +  {y-y'  +  v'Xfy^ ^ | 
xdx  dy  dx'  dy'.  (2.24) 

Equation  (2.24)  can  be  evaluated  numerically  to  give  theoretical  values  for  the  OTF  correla¬ 
tions  between  any  two  arbitrary  spatial  frequencies  {u,v)  and  {u\v').  However,  simulation 
via  random  phase  screen  generator  and  optical  system  models  offers  a  significant  decrease  in 
computational  complexity  over  numerical  evaluation.  A  phase  screen  generator  and  Monte 
Carlo  simulation  are  used  to  generate  the  required  OTF  statistics  in  this  dissertation. 

The  next  section  reviews  classical  estimation  techniques  important  in  processing  turbu¬ 
lence  degraded  images.  These  techniques  include  least  squares  estimation,  Gerchberg-Saxton 
algorithms,  and  ML  estimation.  In  each  case,  basic  theory  and  past  research  related  to  as¬ 
tronomical  imaging  are  reviewed. 
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2.4  Classical  Estimation 

2.4.1  Least  Squares.  Consider  an  inverse  problem  based  on  the  image  model  given 
in  Eqs.  (2.8)  and  (2.9).  A  simple  approach  involves  ignoring  the  noise  and  assuming  no 
a  priori  knowledge  about  the  object.  This  method,  known  as  unconstrained  least  squares 
(ULS),  has  an  objective  function  [43] 

J(6)  =  ||d-J^6||^  (2.25) 


where  d  and  o  are  vector  versions  of  the  detected  image  and  object,  respectively.  The  matrix 
H  is  a.  block-circulant  matrix  representing  the  shift-invariant  PSF.  The  notation  ||  •  |p 
represents  the  Euclidean  norm  of  a  column  vector,  i.e.  ||a||^  =  a^a,  and  the  superscript 
T  denotes  the  matrix  transpose  operator.  The  estimate  that  minimizes  Eq.  (2.25)  is  the 
well-known  expression  [43] 

6=(^H^Hy^H^d  (2.26) 

when  the  matrix  H  is  full  rank,  thus  the  invertibility  of  H^H  is  guaranteed.  Geometrically, 
6  is  the  orthogonal  projection  of  d  onto  the  subspace  spanned  by  the  columns  of  H,  as  shown 
in  Fig.  2.2.  If  is  a  square  matrix  and  has  full  rank,  Eq.  (2.26)  reduces  to  the  intuitive 
form 

6  =  H-^d.  (2.27) 


Recall  that  if  is  a  block-circulant  matrix  and  is  diagonalized  by  the  discrete  Fourier  trans¬ 
form.  Thus,  the  Fourier  domain  equivalent  to  Eq.  (2.27)  is  the  scalar  inverse  filter  which  can 


be  expressed  as  [41] 


d{u,v) 


D{u,v) 

7{{u,v)' 


(2.28) 


Consistent  with  the  original  objective  function  given  at  Eq.  (2.25),  the  inverse  filter  incor¬ 
porates  no  prior  knowledge  about  the  object  or  noise.  In  addition,  Eq.  (2.28)  is  not  valid 
at  spatial  frequencies  {u,v)  where  H{u,v)  =  0.  Even  when  the  OTF  is  non-zero,  the  filter 
amplifies  noise  at  high  spatial  frequencies  [41].  Despite  these  drawbacks,  a  modified  inverse 
filter  is  commonly  applied  to  deconvolve  AO  compensated  images  [72].  An  iterative  imple- 
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(a)  (b) 


Figure  2.2  Geometrical  interpretation  of  ULS  estimation,  (a)  Subspace  spanned  by  the 
linearly  independent  column  vectors  of  the  full  rank  matrix  H.  (b)  The  object 
estimate  6  is  the  orthogonal  projection  of  the  detected  image  d  onto  the  sub¬ 
space  shown  in  (a).  The  quantity  e  denotes  the  error  between  data  and  object 
estimate. 


mentation,  known  as  the  Van  Cittert  method  [2],  is  also  available.  The  iterative  approach 
can  be  advantageous  for  two  reasons:  (l)  the  iterations  can  be  stopped  before  convergence 
and  excessive  noise  amplification,  and  (2)  no  matrix  inversions  are  required  [2] . 

Due  to  the  ill-conditioned  nature  of  the  previous  inversion  problem,  the  ULS  solution 
is  corrupted  by  high  spatial  frequency  noise.  Constraints  can  be  used  to  improve  estimator 
performance  by  incorporating  a  priori  knowledge.  Some  signal  processing  applications  use 
rigid  linear  constraints  which  reduce  the  size  of  the  solution  subspace.  A  typical  linear 
constraint  can  be  written  as  [43] 

Ao  =  b,  (2.29) 

where  A  is  a  known  full  rank  matrix  and  b  is  a  known  vector.  Figure  2.3  illustrates  the 
geometrical  interpretation  associated  with  CLS.  Note  that  the  constrained  solution  6c  can 
be  viewed  as  a  projection  of  the  unconstrained  solution  6  onto  the  constraint  subspace  [43]. 
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Figure  2.3  Geometrical  interpretation  of  CLS  estimation  with  a  rigid  linear  constraint. 

The  constrained  solution  6c  is  a  projection  of  the  unconstrained  solution  6 
onto  the  constraint  subspace. 

Linear  constraints  provide  a  useful  illustration  of  CLS  estimation.  However,  rigid  prior 
knowledge  of  the  true  object  may  not  be  available  when  processing  astronomical  images. 
Instead,  a  flexible  constraint  such  as  smoothness  can  be  used  with  Lagrangian  minimization. 
Hunt  [37]  incorporated  a  quantitative  smoothness  measure  by  minimizing 

IIColP.  (2.30) 

subject  to 

||d-H6|P  =  E{||n|p},  (2.31) 

where  n  is  the  measurement  noise  and  C  is  a  matrix  incorporating  a  smoothness  measure 
such  as  the  two  dimensional  Laplacian  [21].  A  straightforward  Lagrangian  minimization 
yields  the  solution  [37] 

6  =  {h'^H  +  'yC^C)  (2.32) 
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where  7  =  1/A  and  A  is  a  Lagrange  multiplier.  Equation  (2.32)  is  identical  to  Tikhonov- 
Miller  regularization  [59,86].  Thus,  the  CLS  estimator  is  valid  for  a  space- variant  PSF  and 
corresponding  iterative  techniques  similar  to  the  Van  Cittert  method  are  available  [2].  When 
H  and  C  are  block  circulant,  the  Fourier  domain  version  of  Eq.  (2.32)  takes  the  convienent 
scalar  form  [37] 

.  _  H*(u,v)D{u,v) 

“’^^“[7i(«,u)|2  +  7]C(«,t;)r 

where  C(«,  v)  is  the  Fourier  transform  of  the  smoothness  measure.  The  smoothness  measure 

provides  regularization  which  is  controlled  by  7.  The  optimal  Lagrange  multiplier  A  can  be 

found  using  a  Newton-Raphson  minimization  based  on  the  statistics  of  the  noise  process  n 

in  Eq.  (2.31)  [37].  Thus,  a  practical  version  of  the  CLS  filter  in  Eq.  (2.33)  is  iterative.  The 

filter  provides  fidelity  to  the  “rough”  inverse  filter  solution  while  satisfying  the  “smooth” 

constraint. 

There  are  many  variants  on  CLS  estimation  that  have  been  used  in  image  recon¬ 
struction  to  include  weighted  least  squares  [87],  constrained  total  least  squares  [36],  and 
regularized  constrained  total  least  squares  [58].  Recently,  a  new  CLS  estimation  algorithm 
addressed  optimal  use  of  object  model  information  [71].  In  practical  applications,  the  object 
model  could  be  a  low  resolution  image  of  the  object.  This  work,  also  known  as  model-based 
CLS,  is  important  because  it  provides  a  simple  method  for  incorporating  object  model  in¬ 
formation  using  a  CLS-based  algorithm.  CLS  estimation  has  also  been  applied  to  images 
degraded  by  random  blur  [3].  As  of  this  writing,  CLS  estimation  has  not  been  applied  to 
noisy  DWFS  data. 

Least  squares  estimation  provides  a  set  of  straightforward  techniques  which  can  in¬ 
corporate  some  a  priori  knowledge  about  the  imaging  scenario.  In  general,  no  probabilistic 
assumptions  are  made  about  the  data.  In  many  practical  applications,  some  sort  of  iterative 
technique  is  needed  to  find  a  useful  solution.  Another  class  of  iterative  techniques  which  does 
not  require  a  priori  probabilistic  information  about  the  data  will  now  be  discussed.  These 
techniques  are  based  on  the  Gerchberg-Saxton  phase  recovery  algorithm  [19]. 
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2.4-2  Gerchberg- Saxton  Algorithms.  In  1972,  Gerchberg  and  Saxton  introduced  a 
simple  iterative  algorithm  for  the  recovery  of  a  complex  wave  function  from  only  modulus 
measurements  [19].  The  method  depends  on  the  Fourier  transform  relationship  given  by 
the  Van  Cittert-Zernike  Theorem  [23] .  The  input  data  are  the  modulus  measurements  from 
source  and  pupil  planes  as  might  be  available  in  electron  microscopy  [19].  The  algorithm 
begins  with  an  initial  guess  of  the  wave  function  phase.  This  phase  guess  is  combined  with 
the  measured  source  plane  modulus  data  and  then  Fourier  transformed  to  the  pupil  plane 
domain.  Here,  the  measured  pupil  plane  modulus  data  is  imposed  on  the  new  function  and 
the  result  inverse  Fourier  transformed  back  to  the  source  domain.  The  measured  amplitude 
data  is  imposed  again  in  the  source  domain  and  the  process  repeated  until  a  suitable  conver¬ 
gence  criterion  is  met.  The  algorithm  is  based  on  the  fact  that  a  change  in  amplitude  alone 
in  one  domain  of  the  Fourier  transform  results  in  a  change  in  both  amplitude  and  phase 
distributions  in  the  other  domain  [19].  Fienup  [7]  modified  the  Gerchberg-Saxton  algorithm 
to  address  the  phase  retrieval  problem  in  speckle  imaging  [47].  The  problem  here  is  to  find 
an  object  consistent  with  measured  Fourier  modulus  data.  In  finding  an  estimate  of  the 
object,  the  Fourier  domain  phase  is  “retrieved”  in  the  process.  The  Fienup  solution,  known 
as  the  error  reduction  algorithm,  involves  Fourier  transforming  back  and  forth  between  do¬ 
mains,  satisfying  the  constraints  in  one  before  returning  to  the  other.  Figure  2.4  shows  a 
block  diagram  of  this  simple  technique.  The  technique  and  related  extensions,  such  as  the 
input-output  algorithm  [7],  have  been  applied  to  turbulence-degraded  images  [8]. 

The  general  Gerchberg-Saxton  approach  can  be  applied  to  any  problem  in  which  par¬ 
tial  constraints  (data  or  a  priori  information)  are  available  in  each  of  two  domains  [9].  Thus, 
problems  such  as  blind  deconvolution  can  be  addressed  via  this  paradigm.  The  term  blind 
deconvolution  refers  to  an  image  reconstruction  problem  in  which  the  degrading  PSF  is  un¬ 
known.  Figure  2.5  gives  the  block  diagram  of  a  blind  deconvolution  algorithm  first  proposed 
by  Ayers  and  Dainty  [1].  In  this  application,  no  knowledge  is  available  about  the  PSF  or 
object  except  the  degraded  data  and  image  domain  constraints.  Two  modified  inverse  filters 
are  used  to  find  an  object  and  PSF  simultaneously  that  are  consistent  with  the  data  and 
constraints.  A  similar  algorithm  has  been  investigated  which  replaces  the  modified  inverse 
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Initial  Guess 


Figure  2.4  Block  diagram  associated  with  the  error  reduction  algorithm. 

filters  with  pseudo- Wiener  filters  [5] .  This  algorithm  has  also  been  successfully  extended  to 
incorporate  multiple  data  frames  [60] . 

Gerchberg-Saxton  algorithms  provide  a  powerful  alternative  to  least  squares  processing 
in  many  imaging  applications.  However,  as  with  most  iterative  algorithms  involving  non- 
convex  objective  functions,  convergence  to  a  global  “best”  solution  is  not  guaranteed.  The 
error  reduction  algorithm,  as  applied  to  the  phase  retrieval  problem,  is  known  to  stagnate 
after  a  few  iterations  [9].  Solution  uniqueness  is  also  a  concern  [78].  Iterative  blind  decon¬ 
volution  based  on  the  Gerchberg-Saxton  approach  is  also  known  to  suffer  from  convergence 
problems  and  sensitivity  to  the  choice  of  an  initial  estimate  [46].  A  classical  method  which 
introduces  statistical  assumptions  about  the  data  is  now  discussed.  Here,  statistical  analysis 
will  be  important  to  understanding  algorithm  performance. 

2.4-3  Maximum  Likelihood.  Let  us  now  consider  an  image  reconstruction  problem 
in  which  the  probability  density  function  (PDF)  of  the  detected  data  d,  denoted  Pd(d;  o,  H), 
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Figure  2.5  Block  diagram  associated  with  Ayers-Dainty  iterative  blind  deconvolution. 

is  known.  Here,  the  semi-colon  indicates  that  the  PDF  is  parameterized  by  o  and  H.  The 
ML  principle  can  be  stated  as  [43]: 

Given  a  realization  of  d  and  the  PSF  matrix  H,  find  6  which  maximizes  Pd(d;  o,  H). 

In  general,  the  ML  estimate  is  obtained  by  evaluating  the  likelihood  function  using  the 
data  realization  d  and  then  searching  for  the  appropriate  6.  When  the  likelihood  function 
is  continuously  differentiable  in  o  and  convex,  the  ML  estimate  may  be  determined  by 
differentiating  pd(d;o,  ff)  with  respect  to  o,  setting  equal  to  the  zero  vector,  and  solving 
for  o.  ML  estimation  is  a  popular  technique  because  it  is  a  “turn-the-crank”  procedure  [43] 
for  complicated  estimation  problems.  ML  estimation  is  also  closely  related  to  least  squares 
when  the  likelihood  function  is  Gaussian  distributed.  In  this  case,  the  ML  and  weighted 
least  squares  estimators  are  identical  when  the  weighting  matrix  is  chosen  as  the  inverse  of 
the  data  covariance  matrix  [43]. 

In  some  cases,  the  ML  estimator  cannot  be  found  by  taking  the  analytic  derivative  of 
the  likelihood  function.  Another  advantage  of  ML  estimation  is  that  6  can  always  be  found 
by  maximizing  the  likelihood  function  numerically.  A  straightforward  grid  search  can  be  used 
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if  o  is  known  to  exist  on  a  finite  interval.  However,  when  processing  imagery,  the  unknown 
parameters  represent  irradiance  values,  which  do  not  always  exist  on  the  required  finite 
interval.  Thus,  iterative  maximization  methods  must  be  used.  Some  commonly  used  iterative 
maximization  techniques  are  Newton- Raphson  [40],  gradient  descent  [41],  and  expectation- 
maximization  (EM)  [6]  algorithms.  The  EM  algorithm  is  of  particular  interest  in  processing 
astronomical  images  because  it  is  well  suited  to  vector  parameters.  The  EM  algorithm  is 
built  on  the  hypothesis  that  some  data  sets  allow  easier  determination  of  the  ML  estimate 
than  the  data  d.  This  new  data  set  y  is  known  as  the  complete  data,  while  d  is  known  as  the 
incomplete  data.  The  standard  procedure  is  to  assume  that  there  is  a  many-to-one  complete 
to  incomplete  data  transformation.  Thus,  the  EM  strategy  is  to  trade  the  difficult  problem 
of  maximizing  Pd(d;  o,  H)  for  the  easier  problem  of  maximizing  Py(y;  o,  H).  Since  y  does 
not  really  exist,  the  algorithm  incorporates  the  following  iterative  expectation-maximization 
procedure: 

1.  Expectation  (E)  Evaluate  Ey|d[py(y;oj(;,H’)]  using  the  object  estimate.  The  no¬ 

tation  Ea\b  denotes  the  conditional  expected  value  of  the  random  variable  A  given 
B. 

2.  Maximization  (M)  Use  the  conditional  expectation  from  the  previous  step  to 
generate  o^+i. 

3.  Repeat  the  expectation  step  using  o^+i. 

The  EM  algorithm  is  guaranteed  to  converge  to  at  least  a  local  minimum  under  some  mild 
mathematical  conditions  [6,43]. 

ML  estimation  in  all  its  various  forms  has  been  widely  applied  in  image  reconstruction 
to  include  astronomy  [51,68,81],  tomography  [79],  fluorescence  microscopy  [32,34],  and  a 
variety  of  blind  deconvolution  problems  [33,48,77,85].  As  noted  above,  iterative  maximiza¬ 
tion  techniques,  such  as  the  EM  algorithm,  are  widely  used  to  generate  the  ML  estimate. 
These  iterative  techniques  are  not  guaranteed  to  converge  to  a  global  maximum  [43] .  Instead, 
ML  algorithms  may  suffer  convergence  problems  due  to  numerical  inaccuracies,  sensitivity 
to  local  minima,  and  the  choice  of  an  initial  estimate.  Convergence  problems  are  a  special 
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concern  in  the  blind  deconvolution  application  [46].  Finally,  ML  estimation  requires  the 
maximization  of  a  random  likelihood  function.  Thus,  predictive  performance  analysis  can 
be  difficult. 

ML  estimation  is  a  powerful  statistical  technique  which  relies  on  knowledge  of  the 
data  PDF.  In  the  section  below,  Bayesian  estimation,  which  relies  on  statistical  assumptions 
about  data  and  object,  is  discussed.  The  Bayesian  techniques  reviewed  in  this  chapter  include 
the  Wiener  filter,  Kalman  filter,  and  MAP  estimation.  Once  again,  basic  theory  and  past 
research  related  to  astronomical  imaging  are  reviewed. 


2.5  Bayesian  Estimation 

2.5.1  Wiener  Filter.  The  main  drawback  to  ULS  processing  and  the  inverse  filter 
is  sensitivity  to  noise.  This  deficiency  is  no  surprise  since  ULS  does  not  take  into  account 
noise  effects.  In  contrast,  the  Wiener  filter  incorporates  a  statistical  description  of  object 
and  noise  such  that  the  object  estimate  6  can  be  written  as  [64] 

6  =  RooH^iHRooH^  +  i?nn)"'d,  (2.34) 


where  Roo  =  E[oo^]  is  the  object  correlation  matrix  and  Rnn  —  E[nn^]  is  the  noise  corre¬ 
lation  matrix.  In  general  image  reconstruction  applications,  the  mean  object,  denoted  o,  is 
non-zero.  Equation  (2.34)  is  derived  based  on  minimizing  mean  square  error  (MSE)  between 
the  true  object  irradiance  o  and  6.  The  object  estimate  6  is  constrained  to  be  linear  to  the 
related  detected  data  d,  thus  producing  a  mathematically  tractable  non-iterative  solution. 
In  this  application,  the  object  statistical  model  is  static.  A  more  general  form  could  allow 
for  a  dynamic  object  model  [43].  When  object  and  noise  are  stationary,  Rgg  and  Rnn  are 
block  Toeplitz  matrices  [62].  These  correlation  matrices  can  be  made  to  approximate  block 
circulant  matrices  and,  therefore,  are  diagonalized  by  the  discrete  Fourier  transform  [21]. 
The  result  is  a  Fourier  domain  scalar  Wiener  filter  given  as  [21] 


0{u,v) 


H*{u,v)D{u,v) 


(2.35) 
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where  Qn{u^v)  =  E[|A^(it, t;)^]  and  Qo{u,v)  =  E[|0(«, ?;)p]  are  the  noise  and  object  PSDs, 
respectively.  This  approach  has  been  used  in  many  imaging  applications  to  include  enhance¬ 
ment  of  scanning  electron  micrographs  [49],  multichannel  processing  [18],  image  recogni¬ 
tion  [52],  and  deconvolution  of  AO  compensated  images  [72].  The  scalar  Wiener  filter  is 
suboptimal  with  respect  to  MSE  for  non-stationary  object  ensembles.  Since  many  object 
ensembles  are  non-stationary  in  the  mean  and  may  be  non-stationary  with  regard  to  co- 
variance,  the  filter  may  not  incorporate  valuable  a  priori  information  about  a  given  object 
class. 

Pratt  [63]  proposed  a  generalized  vector  Wiener  filter  that  is  optimal  with  respect  to 
MSE  for  non-stationary  object  ensembles  in  signal-independent  noise.  Pratt’s  Wiener  filter 
theory  has  been  extended  to  images  degraded  by  both  blur  and  signal-independent  noise. 
The  resultant  Fourier  domain  filter  expression  can  be  written  as  [76] 

6  =  RooH*^  {Rd.Roo'Hd  +  Rnn}  ^  D,  (2.36) 

where  Rqo  =  E[00^]  is  the  spatial  frequency  correlation  matrix  of  the  object,  R^n  = 
E[NN^]  is  the  spatial  frequency  correlation  matrix  of  the  signal-independent  noise,  Tid 
is  a  diagonal  matrix  of  known  OTF  elements,  and  the  superscript  *  denotes  the  complex 
conjugate  of  the  matrix  elements.  The  role  of  the  object  correlations  between  different  spatial 
frequencies  has  not  been  extensively  studied.  In  addition,  Eq.  (2.36)  does  not  consider  a 
random  OTF  due  to  atmospheric  turbulence. 

Imaging  through  atmospheric  turbulence  is  a  severe  manifestation  of  the  broader  prob¬ 
lem  of  randomness  in  the  pupil  function  of  an  optical  system.  A  random  pupil  function  can 
be  caused  by  a  number  of  factors  such  as  camera  movement  relative  to  an  object,  dust  par¬ 
ticles  on  optical  surfaces,  or  turbulence  in  water.  Wiener  filter  theory  has  been  applied  to 
images  degraded  by  random  blur.  Slepian  [80]  extended  the  Helstrom  scalar  Wiener  filter 
to  incorporate  the  first  order  OTF  statistic  E[H{u,  ?;)]  and  the  second  order  OTF  statistic 
E[|7^(u,  ?;)p]  associated  with  atmospheric  turbulence.  Ward  and  Saleh  [87,88]  investigated 
an  iterative  Wiener  filter  for  one  dimensional  data.  Guan  and  Ward  modified  this  itera¬ 
tive  Wiener  filter  to  process  two  dimensional  images  [25]  and  investigated  a  closely  related 
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geometrical  mean  filter  [26].  This  previous  work  does  not  investigate  the  role  of  OTF  correla¬ 
tions  between  different  spatial  frequencies.  In  addition,  shot  noise  effects  are  not  considered. 
The  Wiener  filters  given  in  Eqs.  (2.34)-(2.36)  incorporate  static  statistical  models.  Both  the 
Wiener  and  Kalman  filters  can  also  incorporate  dynamic  models.  In  the  next  section,  the 
Kalman  filter  in  adaptive  image  reconstruction  is  briefly  discussed. 

2.5.2  Kalman  Filter.  The  Wiener  filter  given  in  Eq.  (2.34)  is  based  totally  on 
knowledge  of  the  data  autocorrelation  and  their  cross-correlation  with  the  object.  Also, 
changes  in  the  object  model  with  respect  to  time  are  not  considered.  General  Wiener- 
Kalman  filter  theory  can  provide  an  adaptive  image  reconstruction  capability  when  dynamic 
statistical  models  are  available.  The  seminal  work  by  Woods  and  Radewan  [96]  led  to  a  two 
dimensional  reduced  update  scalar  Kalman  filter  (RUKF).  Here,  only  the  pixels  in  a  small 
neighborhood  of  the  current  pixel  are  updated.  It  is  assumed  that  a  pixel  is  uncorrelated 
with  other  pixels  outside  this  neighborhood  called  the  update  region  [42].  The  RUKF  has 
been  applied  to  deconvolution- type  problems  [95].  More  recent  advances  in  Kalman  filter 
processing  of  two  dimensional  data  involve  a  reduced  order  model  (ROM)  representation  [42] 
and  the  ROM  Kalman  filter  [97] .  A  two  dimensional  Kalman  filter  has  also  been  investigated 
for  images  degraded  by  random  blur  [67].  However,  the  necessary  dynamic  object  state  model 
is  not  available  in  many  applications. 

The  next  section  introduces  a  Bayesian  technique  which  is  closely  related  to  the  ML 
estimator.  However,  the  MAP  estimator  incorporates  a  prior  PDF  associated  with  the  object. 

2.5.3  Maximum  A  Posteriori.  The  MAP  estimation  principle  can  be  stated  as 
follows  [43]: 

Given  a  realization  of  d  and  the  PSF  matrix  H,  find  6  which  maximizes  the 

posterior  PDF  po(o[d;  H). 

The  posterior  PDF  is  a  conditional  distribution.  Maximizing  Po(o[d;  H)  has  been  shown 
to  minimize  a  “hit-or-miss”  cost  function  which  assigns  no  penalty  for  small  errors  and 
maximum  penalty  for  all  errors  in  excess  of  a  threshold  S  [43].  This  cost  function  can  be 
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expressed  mathematically  as  [43] 


Cost(e) 


1 


<5 

>S 


(2.37) 


where  the  variable  e  represents  error  and  5  >  0.  If  is  made  very  small,  this  cost  function 
assigns  the  same  penalty  for  all  errors  and  no  penalty  for  no  error. 

The  posterior  PDF  can  be  written  in  a  more  intuitive  form.  Using  Bayes  rule  [62], 
Po(o|d;  H)  becomes 


Po{o\d-,H) 


Pd(d;  H\o)po{o) 


(2.38) 


Pd(d;  o,  H) 

Thus,  maximizing  the  numerator  of  Eq.  (2.38)  is  equivalent  to  maximizing  the  posterior 
PDF.  Notice  that  Pd(d;  JT|o)po(o)  is  very  similar  to  the  ML  likelihood  function  given  in 
Section  2.4.3  except  that  the  object  is  now  considered  a  random  process  with  known  prior 
PDF  Po{o).  As  in  the  ML  case,  a  candidate  MAP  estimate  can  be  found  by  evaluating 
Pd(d;  H\o)po{o)  at  the  given  data  realization,  differentiating  with  respect  to  o,  setting  equal 
to  the  zero  vector,  and  solving  for  o.  When  a  closed  form  solution  is  not  practical,  the 
MAP  estimate  can  be  found  via  iterative  maximization  methods  such  as  gradient  ascent 
algorithms  [24,38]  and  the  EM  algorithm  [28].  Iterative  MAP  estimation  can  suffer  from  the 
same  numerical  convergence  problems  noted  for  ML  estimation  in  Section  2.4.3. 


2. 6  Summary 

This  chapter  provided  background  material  related  to  image  reconstruction.  The  em¬ 
phasis  here  is  on  the  deconvolution  of  astronomical  images.  Thus,  the  image  degradation 
model  associated  with  atmospheric  turbulence,  shot  noise,  and  detector  read  noise  was  re¬ 
viewed  as  well  as  statistical  theory  related  to  atmospheric  turbulence.  The  rest  of  the  chapter 
presented  important  classical  and  Bayesian  estimation  techniques,  most  of  which  involve  it¬ 
erative  optimization.  While  these  techniques  are  powerful,  drawbacks  related  to  convergence 
and  ease  of  analysis  do  exist.  Current  linear  reconstruction  methods,  such  as  the  scalar 
Wiener  filter,  do  not  improve  Fourier  domain  signal-to-noise  ratio  (SNR).  In  addition,  the 
application  of  object,  OTF,  and  noise  correlations  between  different  spatial  frequencies  has 
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not  been  investigated.  In  the  rest  of  this  dissertation,  two  techniques  which  bridge  the  gap 
between  linear  processing  and  nonlinear  iterative  optimization  are  presented:  CLS  process¬ 
ing  of  DWFS  data  and  the  vector  Wiener  filter.  The  vector  Wiener  filter  provides  a  useful, 
non-iterative  complement  to  nonlinear  iterative  optimization.  The  complete  vector  Wiener 
filter  is  derived  in  Chapter  IV  with  performance  data  given  in  Chapters  V  and  VI.  First, 
Chapter  III  presents  a  novel  application  of  CLS  post-detection  image  processing. 
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III.  Constrained  Least  Squares  Incorporating  Wavefront  Sensing 

3.1  Introduction 

In  this  chapter,  a  technique  for  processing  noisy  deconvolution  from  wavefront  sens¬ 
ing  (DWFS)  data  based  on  constrained  least  squares  (CLS)  estimation  is  presented.  The 
new  algorithm  selects  a  value  for  the  regularization  constant  which  is  consistent  with  the 
ensemble-averaged  data  and  a  constraint  [11].  This  problem  is  solved  using  the  Lagrange 
multiplier  technique  [21,37].  A  closed  form  solution  for  the  object  estimate  is  obtained  which 
is  analogous  to  the  traditional  DWFS  estimator  [20, 75]  with  the  Lagrange  multiplier  serving 
as  a  regularization  constant.  An  iterative  approach  based  on  a  Newton- Raphson  minimiza¬ 
tion  [40, 71]  is  used  to  find  the  optimal  Lagrange  multiplier.  The  iteration  incorporates  the 
statistics  of  both  the  optical  transfer  function  (OTF)  and  noise.  The  sample  results  given 
here  demonstrate  that  CLS  estimation  provides  high  quality  processing  of  noisy  DWFS  data 
automatically.  No  ad  hoc  tuning  of  the  regularization  constant  is  necessary.  CLS  object 
estimates  are  compared  with  those  processed  via  manual  parameter  selection.  In  all  cases, 
the  new  technique  provides  images  with  comparable  resolution.  In  addition,  the  algorithm 
is  computationally  inexpensive,  converging  to  a  solution  in  less  than  10  iterations  [11]. 

The  rest  of  this  chapter  is  organized  as  follows.  Section  3.2  reviews  the  traditional 
DWFS  estimator.  In  Section  3.3,  a  CLS  algorithm  is  derived  which  incorporates  a  new 
constraint  based  on  the  ensemble-averaged  DWFS  data.  Section  3.4  gives  sample  results 
associated  with  various  imaging  conditions,  while  Section  3.5  provides  a  brief  summary. 

3.2  Traditional  Estimator 

Raw  DWFS  data  consists  of  an  ensemble  of  short  exposure  images  and  corresponding 
estimates  of  the  wavefront  phase  from  a  wavefront  sensor  (WFS).  The  phase  estimates  ^i{x,  y) 
can  be  used  to  generate  an  estimate  of  the  OTF  'Hi{u,v)  via  a  normalized  autocorrelation 
of  the  generalized  pupil  function  as  shown  in  Eq.  (2.17).  The  subscript  i  refers  to  the  i^^ 
realization,  while  the  tilde  denotes  a  quantity  estimated  directly  from  information  provided 
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by  WFS  hardware.  The  DWFS  estimator  suggested  by  Primot  et  al.  [66]  can  be  written  as 


(H*{u,v)I{u,v)) 
0{u,v)  =  - - 


2\  > 


where  {X{u,v))  is  defined  as 


{X{u,v))  =  2'^Xi(u,v), 


i=l 


(3.1) 


(3.2) 


I{u,  v)  is  the  noiseless  image  Fourier  spectrum  and  L  is  the  number  of  images  in  the  ensemble. 
Equation  (3.1)  does  not  provide  acceptable  reconstructions  due  to  residual  noise  in  the  OTF 
estimation  process.  When  detector  noise  is  present,  I(u,v)  is  not  available  and  the  detected 
image  Fourier  spectrum  D(u,v)  must  be  substituted  in  Eq.  (3.1).  The  standard  solution 
is  to  incorporate  a  regularization  constant  e  in  the  filter  denominator  which  gives  the  new 
estimator  [20,93] 


0{u,v) 


(H*{u,v)D{u,v)'^ 


H{u,  v) 


+  6 


(3.3) 


The  regularization  constant  e  is  adjusted  in  a  completely  ad  hoc  manner  by  the  user  based 
on  the  perceived  image  quality. 


In  the  next  section,  a  CLS  algorithm  is  derived  which  takes  advantage  of  noise  reduction 
through  ensemble  averaging  by  directly  incorporating  {'H*{u,v)D{u,v))  and  {\'H{u,v)\^). 
Thus,  the  objective  function  used  to  derive  this  estimator  takes  on  an  unfamiliar  form  when 
compared  to  more  traditional  CLS  applications  in  image  processing  [21,37,71]. 


3.3  Modified  Constrained  Least  Squares  Formulation 

3.3.1  Problem  Statement.  This  CLS  estimation  problem  can  be  stated  as  follows: 
“Given  the  DWFS-derived  estimate  {'H*(u,v)D{u,v)),  the  ensemble-averaged  magnitude- 
squared  OTF  estimate  {\'H{u,v)\'^),  and  statistical  models  for  the  OTF  and  noise,  find  the 
CLS  optimal  estimate  of  the  object  that  caused  the  detected  image  ensemble.”  To  accomplish 
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this  task,  consider  finding  an  object  estimate  6  that  minimizes 


(ll^caip), 


(3.4) 


subject  to  the  constraint 


\\{H'^d)-{H^H)6\\‘^^E{\\H^n\\%  (3.5) 

where  (7  is  a  block-circulant  constraint  matrix  which  enforces  prior  knowledge  of  the  true 
object  such  as  smoothness  or  support  and  Hi  is  the  block-circulant  estimated  PSF  matrix  for 
the  realization.  Here,  the  {'H*(u,v)D(u,v))  and  {\i-C{u,v)\^)  quantities  are  incorporated 
in  the  constraint  function  directly.  Also,  the  MSE  associated  with  the  measurement  is  now 
E{||f/'^n|p}  instead  of  the  norm-squared  of  the  noise  E{||n|p}  used  previously  [21,37].  Even 
though  Eqs.  (3.4)  and  (3.5)  are  unfamiliar  with  respect  to  traditional  CLS  applications  in 
image  processing  [21,37,71],  it  will  be  shown  that  the  associated  Fourier  domain  filter  has 
the  same  general  form  as  Eq.  (3.3). 

3.3.2  Closed  Form.  Solution.  The  appropriate  objective  function  J(6,  A)  for  the 
Lagrange  minimization  can  be  written  as 

J(6,  A)  =  d^C^{H^H)Cd  +  x{\\{H^d)  -  {H^H)6\\'^  -  E{||i7^n||2}}  ,  (3.6) 

where  A  is  a  Lagrange  multiplier  [21,37].  Since  the  Lagrangian  in  Eq.  (3.6)  is  quadratic,  it 
can  be  minimized  by  differentiating  J(6,  A)  with  respect  to  6,  setting  the  derivative  equal  to 
zero,  and  solving  for  6.  Applying  the  appropriate  vector-matrix  identities  [56]  to  take  the 
derivatives  in  Eq.  (3.6)  yields 

6  =  [{H^H){H^H)  +  'yC'^{H^H)Cy'  {H^H){H^d),  (3.7) 

where  7  =  1/A.  It  is  not  computationally  efficient  to  evaluate  Eq.  (3.7)  directly  since  C 
and  H  are  large  matrices  for  standard-sized  image  arrays.  However,  these  matrices  are 
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block  circulant.  Block  circulant  matrices  are  diagonalized  by  the  discrete  Fourier  transform 
allowing  the  transform  domain  equivalent  of  Eq.  (3.7)  to  provide  the  simpler  scalar  form 
[21,37] 


0(u,  v)  = 


(7l*{u,v)D(u,v)'^ 


(3.8) 


where  C{u,v)  is  the  Fourier  transform  of  the  constraint  c{x,y)  associated  with  the  matrix 
C.  When  C{u,v)  =  1  for  all  spatial  frequencies  (u,v),  Eq.  (3.8)  is  identical  to  Eq.  (3.3). 
However,  Eq.  (3.8)  will  accommodate  a  general  Fourier  domain  constraint  function  which 
can  be  tailored  to  prior  knowledge  or  a  specific  application. 


3.S.S  Newton-Raphson  Iterative  Solution.  Instead  of  tuning  7  manually,  a  Newton- 
Raphson  technique  is  used  to  find  the  parameter  value  which  minimizes  the  objective  function 
J(6,  A)  given  in  Eq.  (3.6).  To  implement  the  Newton-Raphson  technique,  the  derivative  of 
J(6,  A)  with  respect  to  A  is  required.  This  derivative  is  straightforward  to  derive  and  given 
by 

=  IKtf^d)  -  {6^6)6^  -  E{||ff^n|p}.  (3.9) 

This  iteration  can  also  be  accomplished  with  respect  to  the  7  parameter  [37].  However,  the 
required  derivative  is  not  as  straightforward  as  Eq.  (3.9)  [37,71].  The  first  norm-squared 
term  in  Eq.  (3.9)  is  a  function  of  the  DWFS  data  and  the  object  estimate.  It  can  be  com¬ 
puted  conveniently  using  the  Fourier  domain  quantities  {'H*{u,v)D{u,v))  and  {\'H{u,v)\'^). 
However,  the  second  norm-squared  term  E{||.H’^n|]^}  is  not  a  function  of  the  DWFS  data 
and  must  be  derived. 

Rewriting  E{||fZ'^n]]^}  in  the  Fourier  domain  using  Parseval’s  Theorem  [21]  yields 


Substituting  the  image  model  given  in  Eq.  (2.7)  into  Eq.  (3.10)  introduces  the  noise  statistics 
such  that 

I  U,V 

K 

-  exp[-j27r(Ma;„  +  vyn)] 

n=l 

9  ^ 

P  ^ 

-  ^  np  exp  {-j27r(ua;p  +  ?;z/p)})  , 

p-i 

/ 

where  {Xn,yn)  is  the  location  of  the  photoevent  in  the  image  plane,  {xp,yp)  is  the  loca¬ 
tion  of  the  image  pixel,  and  K  is  the  average  number  of  photoevents  per  image.  The 
normalized  object  spectrum  0„(«,  u)  =  0(u,v)/K  has  also  been  introduced  to  the  above 
expression.  Now  the  linearity  property  can  be  used  to  move  the  expectation  inside  the  sum¬ 
mation  over  spatial  frequencies  {u,v)  in  order  to  evaluate  one  term  of  Eq.  (3.11).  At  this 
point,  the  right  side  of  Eq.  (3.11)  can  be  simplified  using  standard  techniques  for  evaluat¬ 
ing  expectations  of  doubly  stochastic  Poisson  random  processes  [74].  These  techniques  use 
nested  conditional  expectations  over  the  random  quantities  (xn,2/n)^  K,  and  H  [74].  Nine 
sub-terms  of  Eq.  (3.11)  are  evaluated  below: 

1.  E{{Ky\On{u,v)\‘^\H{u,v)\^}.  Here,  the  only  random  quantity  is  \'H{u,v)\^.  Thus,  the 
K  and  Oniu,  v)  terms  are  brought  outside  the  expectation  yielding  the  final  result 

Term  1  =  {KY\On{u,v)\^E{\'H{u,v)\^].  (3.12) 

2.  -E  {KOn{u,  v) \H{u,  v)\^H{u,  v)  En=i  exp[j27r(ita:„ -h  upn)]}-  The  random  quantities 
underlying  this  expectation  are  (a;„,y„),  K,  and  ?{.  Bayes  rule  [62]  can  be  used  to 
rewrite  the  joint  probability  density  function  (PDF)  in  the  first  term  using  conditional 
PDFs.  These  conditional  PDFs  translate  to  conditional  expectations  such  that  the 
second  term  can  be  written  as 


(3.11) 
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Term  2  =  —KOn{u,v) 


X  {  [Hiu,  v)\‘^H{u,  v)Ek\H  ^  '^Xn,yn\K,H 


r  K 


exp[j27r(Ma;„  +  vy^)] 

In—l 


(3.13) 


The  innermost  expectation  ^xnyyn\K,H[*]  has  been  evaluated  previously  [74].  The  result 
is  repeated  below  as 

=  KH'(u,  «)0;(»,  v),  (3.14) 

where  the  •  notation  represents  the  argument  of  the  innermost  expectation  in  Eq.  (3.13). 
Substituting  Eq.  (3.14)  into  Eq.  (3.13)  and  evaluating  the  remaining  nested  conditional 
expectations  gives  the  final  result  for  the  second  term 

Term  2  =  -(Kf\On{u,v)\^E{\H{u,v)\^}.  (3.15) 


3.  -E  {KOn{u,  v)\H{u,  v)\^H{u,  v)  %  eyi-p\j2'K{uXp  +  t^2/p)]}.  In  Chapter  II,  the  de¬ 
tector  read  noise  random  variable  Up  was  assumed  to  be  zero  mean  and  statistically 
independent  of  K  and  H.  Thus,  the  third  term  is  zero. 

Term  3  =  0.  (3.16) 


4.  -E  [koi(u,  v)\n{u,v) !%•(«,  v)  E*  1  exp[— j27r(MX„  -f  ^^^/n)]}•  Term  4  is  the  complex 
conjugate  of  Term  2  above.  Thus,  the  final  result  for  the  fourth  term  is 

Term  4  =  -{Ky\On{u,v)\^E{\'H{u,v)\^}.  (3-17) 

5.  E  {|7i(M,  v)\^  En=i  Em=i  exp[-j2Tr{u{xn  -  Xm)  +  v{yn  -  ym))]}-  The  random  quanti¬ 
ties  underlying  this  expectation  are  {xn,yn),  {xmiym),  K,  and  K.  As  with  the  second 


3-6 


term  above,  conditional  expectations  can  be  used  to  rewrite  the  fifth  term  such  that 


Term  5  =  £•«  {^H{u,v)\^'Ek\h  {^xn,yn,xm,ym\K,H  [•]})  •  (3-18) 

Here,  •  represents  the  double  summation  shown  in  the  first  line  of  5  above.  The 
innermost  expectation  ^xn,vn,xm.,ym\K,H[*\  been  evaluated  previously  [74].  The  result 
is  repeated  below  as 


+  (K^  -  K)\On{u,v)\^\n{u,v)\\  (3.19) 

To  evaluate  the  expectation  over  K,  recall  that  the  image  photon  count  is  Poisson 
distributed  which  implies  that  =  {Ky  +  K  [62,74].  Thus,  the  final  result  for  Term 
5  is 

Term  5  =  KE{\n{u,v)\‘^}  +  (Ky\On{u,v)\^E{\n{u,v)\^}.  (3.20) 

6.  E  {|'H(«,  v)\^  (En=i  exp[-j27r(«rc„  +  uy„)])  (Ep=i  np  exp[j27r(uxp  + 1;?/^)])  }.  The  sixth 
term  is  zero  since  Up  is  zero  mean  and  statistically  independent  of  {xn,yn)j  K,  and  7i. 

Term  6  =  0.  (3.21) 

7.  -E  {FO;(«,  v)\H{u,  v)\^H*{u,  v)  E^=i  rip  exp  [—j27r(ua;p  +  vj/p)]}-  The  seventh  term  is 
zero  since  rip  is  zero  mean  and  statistically  independent  of  K  and  7Y. 

Term  7  =  0.  (3.22) 

8.  E{|W(«,»)p(Ej=.n,  exp[-j27r(ita;p  +  uyp)])  exp[j27r(ua;n  +  ^^Z/n)])}-  The  eighth 

term  is  the  complex  conjugate  of  the  sixth  term.  Thus, 

Term  8  =  0.  (3.23) 
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9.  E^\H{u,  v)\^  Ep=i  T,g=i  npUq  exp[-j27r{u(xp  -  Xg)  +  v{yp  -  2/g))]}.  First,  recall  that  H 
and  Up  are  statistically  independent.  Therefore,  the  ninth  term  can  be  written  as 


Term  9  =  E  { \H{u,  'y)  T}  E  |  5^  ripUg  exp[-j2Tv{u{xp  -  Xg)  +  v{yp  -  y,))]  |  (3.24) 


q-l  ) 

The  second  expectation  in  Eq.  (3.21)  has  been  evaluated  previously  [74]  for  the  as¬ 
sumptions  given  in  Chapter  II  such  that  the  final  result  for  Term  9  can  be  expressed 
as 


Term  9  =  P(tIE{\'H{u,v)\‘^}. 


(3.25) 


Adding  the  nine  sub-terms  presented  above  gives  the  result  for  a  single  spatial  frequency 
ofEq.  (3.11): 

E{\n*{u,v)N{u,v)\^}  =  (K  +  Pal)E{\Hiu,v)\‘^}.  (3.26) 

Summing  the  result  in  Eq.  (3.26)  over  all  spatial  frequencies  gives  the  final  result 

E{||F’’n|n  =  +  y;E{|«(»,»)p}  .  (3.27) 

UiV 

If  the  OTF  statistics  are  unavailable  in  Eq.  (3.27),  the  OTF  estimate  data  can  be  substituted, 
which  yields 

E{||ff^n|p}»  (F+Pa,*)X:{|H(«,.,)p).  (3.28) 

UjV 

The  result  given  in  Eq.  (3.28)  is  used  in  Eq.  (3.9)  to  find  the  current  derivative  of  J(6,  A)  with 
respect  to  the  Lagrange  multiplier  A.  This  derivative  is  then  used  in  the  Newton-Raphson 
iteration  to  update  the  value  of  A  and  generate  an  object  estimate  using  Eq.  (3.7)  or  its 
Fourier  domain  equivalent  Eq.  (3.8).  The  iteration  continues  until  the  object  estimate  meets 
a  pre-determined  stopping  criterion.  A  variety  of  criteria  can  be  used  to  stop  the  Newton- 
Raphson  iteration  [40].  Since  the  goal  is  to  minimize  the  Lagrangian  given  in  Eq.  (3.6),  the 
iteration  is  stopped  when  J(6,  A)  is  sufficiently  small.  All  CLS  images  shown  in  this  chapter 
were  processed  until  J(6,  A)  <  0.001.  For  applications  in  which  the  appropriate  value  for  the 
Lagrangian  is  unknown,  the  algorithm  can  terminate  when  the  change  in  A  has  stabilized 
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from  iteration  to  iteration.  In  the  next  section,  sample  results  are  presented  to  illustrate  the 
technique. 

3.4  Sample  Results 

In  this  section,  sample  results  are  presented  which  illustrate  CLS  processing  of  DWFS 
data.  A  variety  of  imaging  conditions  are  investigated  by  varying  turbulence  strength,  light 
level,  and  detector  read  noise  variance.  Data  is  presented  for  uncompensated  and  adaptive 
optics  (AO)  compensated  images  of  a  representative  satellite  object. 

3.4-1  Assumptions.  The  satellite  object  shown  in  Fig.  3.1  was  used  to  generate  all 
data  in  this  chapter.  The  detector  array  was  256  x  256  elements  for  a  total  of  P  =  65536 
pixels.  The  satellite  was  10m  in  length  and  in  low  earth  orbit  at  a  range  of  500km.  The 
diameter  of  the  telescope  was  Im  with  both  imaging  and  wavefront  sensor  wavelengths  set 
at  A  =  500nm.  The  simulated  spectral  bandwidth  was  ±5  %  of  A,  with  the  object  assumed 
to  have  the  same  spectral  signature  as  the  sun.  The  WFS  subaperture  size  was  10cm  for  a 
total  of  60  subapertures  within  the  entrance  pupil.  The  AO  system  model  had  1.2  actuators 
per  To. 

Atmospheric  turbulence  and  detector  noise  effects  were  modeled  using  an  existing  com¬ 
puter  simulation  [74].  For  each  data  realization,  the  simulation  created  a  random  phase 
screen  with  the  appropriate  turbulence  statistics  [89],  calculated  the  true  OTF,  and  formed 
a  detected  image  realization  di{x,y).  At  the  same  time,  a  WFS  model  was  used  to  generate 
a  phase  estimate  ^i{x,y),  which  was  then  used  to  compute  the  estimated  OTF.  Finally, 
the  required  quantities  ii*(u,v)Di{u,v)  and  \'Hi{u,v)\'^  were  accumulated  and  the  process 
repeated  150  times  to  generate  the  ensemble  averages.  The  simulation  incorporates  an  in¬ 
tensity  splitter  which  sends  40%  of  the  photons  to  the  image  plane  and  60%  to  the  WFS. 
Integration  times  of  10ms  were  used  for  both  WFS  and  imaging  system.  A  range  of  object 
brightness  was  modeled  using  four  visual  magnitude  levels  =  +2,  +4,  -|-6,  and  -|-8.  The 
visual  magnitude  allows  astronomers  to  compare  object  brightness  in  the  night  sky  [57]. 
A  step  in  visual  magnitude  indicates  a  corresponding  factor  2.5  change  in  brightness,  with 
smaller  values  of  mj,  indicating  brighter  objects  [57].  The  resulting  photoevents  per  inte- 
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Figure  3.1  Satellite  computer  rendering  used  to  test  CLS  algorithm  performance.  Negative 
image  shown  for  clarity. 

Table  3.1  Average  Photoevents  per  Integration  Time. 


gration  time  for  these  cases  are  presented  in  Table  3.1,  where  ATvv  is  the  average  number  of 
photoevents  across  an  individual  WFS  subaperture. 

3.4-2  Uncompensated  Images.  Figure  3.2  gives  the  detected  short  exposure  im¬ 
age  and  CLS  algorithm  output  associated  with  excellent  seeing  conditions  {tq  =  20cm),  a 
moderate  light  level  satellite  object  {m^  =  +4),  and  detector  read  noise  representing  a  high 
quality  CCD  detector  {a^  =  15  electrons  per  pixel).  A  single  short  exposure  image  is  given 
in  (a)  to  illustrate  shot  noise  and  detector  read  noise  effects.  Clearly,  this  image  is  degraded 
by  the  detection  process.  After  CLS  processing  of  the  DWFS  data  as  shown  in  (b),  image 
resolution  is  greatly  enhanced.  In  this  case,  the  CLS  algorithm  automatically  selected  the 
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Figure  3.2  CLS  algorithm  comparison  with  manual  parameter  selection,  Tq  =  20cm,  = 

+4,  (Tr  =  15  electrons  per  pixel,  (a)  Single  short  exposure  image,  (b)  CLS 
algorithm  estimate,  7  =  0.0015,  4  iterations,  (c),  (d),  (e),  and  (f)  manual 
parameter  selection,  e  =  0.1,0.01,0.0001,  and  0.00001,  respectively. 
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regularization  parameter  7  =  0.0015  in  4  iterations.  In  (c),  (d),  (e),  and  (f)  the  regular¬ 
ization  parameter  e  associated  with  the  traditional  DWFS  estimator  given  in  Eq.  (3.3)  was 
selected  manually.  For  images  (c)  and  (d),  e  >  7.  For  images  (e)  and  (f),  e  <  7.  Notice 
that  the  quality  of  the  CLS  image  in  (b)  is  as  good  or  better  than  the  manual  images.  Thus, 
Fig.  3.2  illustrates  that  the  CLS  algorithm  selected  a  reasonable  value  for  7  in  this  case. 

Now  consider  a  brighter  object  and  different  seeing  conditions.  Figure  3.3  gives  the 
detected  image  and  CLS  algorithm  output  associated  with  changing  turbulence  strength  and 
a  brighter  satellite  object  (mj,  =  -t-2).  The  CCD  detector  read  noise  remains  unchanged  from 
the  previous  case.  Images  (a)  and  (c)  provide  the  detected  image  data  for  the  To  =  10cm 
and  20cm  cases,  respectively.  Images  (b)  and  (d)  give  the  corresponding  CLS  algorithm 
output.  Here,  7  :=  0.000012  for  the  To  =  10cm  case  and  7  =  0.000067  for  the  =  20cm 
case.  The  CLS  algorithm  provides  more  regularization  or  “smoothing”  as  To  increases.  This 
observation  is  consistent  with  the  form  of  Eq.  (3.28),  where  better  seeing  conditions  lead  to 
larger  values  for  the  quantity  (|7t((m,u)P^.  In  general,  larger  values  for  the  quantity 
Y,u,v  (|^(W)'*^)P)  Isad  to  larger  7  values  and  more  noise  smoothing. 

While  fo  influences  CLS  algorithm  performance,  object  brightness  provides  a  more 
severe  limit.  Not  only  does  shot  noise  degrade  the  short  exposure  data,  it  also  restricts  WFS 
accuracy.  Without  a  sufficiently  accurate  OTF  estimate,  DWFS  performance  is  severely 
degraded.  To  illustrate  these  limitations,  consider  Fig.  3.4.  Here,  To  =  10cm  and  detector 
read  noise  remains  unchanged  from  the  previous  cases.  Images  (a),  (c),  and  (e)  show  the 
short  exposure  data  for  object  brightness  cases  m„  =  +4,  -|-6,  and  +8,  respectively.  In  (c) 
and  (e),  noise  dominates  the  data  realization  to  the  point  that  no  satellite  image  is  visible. 
Images  (b),  (d),  and  (f)  give  the  corresponding  CLS  estimates  where  7  =  0.00027,0.0093, 
and  0.8180,  respectively.  As  object  brightness  decreases,  the  output  images  are  more  blurred. 
This  effect  is  consistent  with  a  relatively  large  7  value  and  deconvolution  using  a  poor  quality 
estimate  of  the  OTF. 

To  emphasize  the  limitations  imposed  by  shot  noise  further,  consider  Fig.  3.5.  Here, 
ro  =  10cm  and  =  -1-4.  Detector  read  noise  variance  is  adjusted  such  that  image  (a)  was 
collected  with  a  low  noise  array  (<7^  =  10  electrons  per  pixel)  and  image  (b)  with  a  high 
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Figure  3.3  CLS  algorithm  output  versus  atmospheric  turbulence  strength,  rrii,  =  +2,  ar  = 
15  electrons  per  pixel,  (a)  Detected  image,  ro  =  10cm.  (b)  CLS  estimate. 
To  =  10cm,  7  =  0.000012,  4  iterations,  (c)  Detected  image.  To  =  20cm.  (d) 
CLS  estimate,  Vg  =  20cm,  7  =  0.000067,  4  iterations. 


Figure  3.4  CLS  algorithm  output  versus  object  brightness,  Vg  =  10cm,  ar  =  15  electrons 
per  pixel,  (a)  Detected  image,  =  +4.  (b)  CLS  estimate,  m,,  =  +4,  7  = 
0.00027,  4  iterations,  (c)  Detected  image,  =  +6.  (d)  CLS  estimate,  = 
+6,  7  =  0.0093,  5  iterations,  (e)  Detected  image,  m„  =  +8.  (d)  CLS  estimate, 
m„  =  +8,  7  =  0.8180,  7  iterations. 
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Figure  3.5  CLS  algorithm  output  versus  detector  read  noise  strength,  Vg  =  10cm,  rrii,  = 
+4.  (a)  Detected  image,  =  10  electrons  per  pixel,  (b)  CLS  estimate,  =  10 
electrons  per  pixel  ,  7  =  0.00014,  4  iterations,  (c)  Detected  image,  cr^  =  30 
electrons  per  pixel,  (d)  CLS  estimate,  cr^  =  30  electrons  per  pixel ,  7  =  0.00097, 
4  iterations. 

noise  array  {cr  =  30  electrons  per  pixel).  Images  (b)  and  (d)  give  the  corresponding  CLS 
estimates  where  7  =  0.00014  and  0.00097,  respectively.  Notice  that  an  increase  in  read  noise 
strength  does  not  have  the  drastic  effect  on  image  quality  observed  in  Fig.  3.4.  Detector 
read  noise  does  not  affect  the  accuracy  of  the  OTF  estimate  provided  by  the  WFS.  Thus, 
the  effect  on  algorithm  output  is  less  pronounced. 

3.4-3  Adaptive  Optics  Compensated  Images.  DWFS  processing  can  also  be  applied 
to  AO  compensated  images.  Here,  the  algorithm  deconvolves  atmospheric  turbulence  effects 
associated  with  the  residual  error  between  the  true  phase  perturbation  and  the  phase  imposed 
by  the  deformable  mirror  (DM).  As  noted  earlier,  the  ratio  between  DM  actuator  number 
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(c)  (d) 

Figure  3.6  CLS  algorithm  output  for  AO  compensated  images,  rrii,  =  +4,  cr^  =  15  electrons 
per  pixel,  (a)  Detected  image,  Tg  —  10cm.  (b)  CLS  estimate,  Tg  =  10cm, 
7  =  0.0048,  4  iterations,  (c)  Detected  image,  Vg  =  20cm.  (d)  CLS  estimate, 
Tg  =  20cm,  7  =  0.0039,  4  iterations. 
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and  To  was  fixed  at  1.2.  Thus,  the  number  of  actuators  will  change  as  To  is  adjusted.  Two 
AO  systems  associated  with  Tq  =  10cm  (57  DM  actuators)  and  Tq  =  20cm  (13  DM  actuators) 
are  modeled.  Figure  3.6  illustrates  CLS  algorithm  performance  on  AO  compensated  images 
when  rUi,  =  +4  and  detector  read  noise  is  =  15  electrons  per  pixel.  Images  (a)  and  (c) 
give  examples  of  raw  AO  compensated  data  when  Tq  =  10cm  and  20cm,  respectively.  Images 
(b)  and  (d)  show  the  corresponding  CLS  estimates,  where  7  =  0.0048  in  (b)  and  7  =  0.0039 
in  (d).  The  CLS  algorithm  does  a  good  job  of  noise  suppression  for  AO  compensated  images 
for  the  same  reason  noted  above  in  the  discussion  associated  with  Fig.  3.3.  Here,  AO 
compensation  leads  to  larger  values  for  the  quantity  J2u,v  Eq.  (3.28)  and 

larger  7  values. 

3.5  Summary 

A  CLS  estimator  that  incorporates  noisy  DWFS  data,  noise  statistics,  and  OTF  statis¬ 
tics  was  investigated.  For  a  particular  choice  of  Fourier  domain  constraint,  the  estimator 
selects  the  regularization  parameter  automatically.  No  ad  hoc  tuning  is  necessary  to  reduce 
high  spatial  frequency  noise  effects  in  the  DWFS  image.  The  CLS  estimator  uses  a  Newton- 
Raphson  iteration  to  select  a  Lagrange  multiplier  which  minimizes  an  objective  function. 
The  objective  function  uses  ensemble-averaged  data  directly,  which  aids  in  noise  suppres¬ 
sion.  The  sample  results  show  that  the  new  algorithm  produces  DWFS  images  comparable 
in  quality  to  manual  regularization  with  minimal  computational  expense.  While  turbulence 
and  detector  read  noise  strength  impact  algorithm  performance,  shot  noise  imposes  a  fun¬ 
damental  limit  on  the  deconvolution  process.  Finally,  the  CLS  estimator  was  derived  to 
incorporate  a  general  Fourier  domain  constraint.  Thus,  other  constraint  functions  can  be 
used  based  on  the  specific  application. 

In  the  next  chapter,  a  non-iterative  Bayesian  deconvolution  filter  is  derived.  Unlike 
this  CLS  approach,  the  vector  Wiener  filter  incorporates  object,  OTF,  and  noise  correlations 
between  different  spatial  frequencies.  Chapters  V  and  VI  will  further  explore  the  value  of 
this  model-based  information  for  several  applications  related  to  astronomical  imaging. 
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IV.  Vector  Wiener  Filter  Analysis 


4.1  Introduction 

In  this  chapter,  a  Fourier  domain  vector  Wiener  filter  is  derived  which  incorporates 
complete  object,  blur,  and  noise  correlation  statistics.  The  derivation  extends  the  original 
work  by  Pratt  [63]  to  account  properly  for  a  random  optical  transfer  function  (OTF)  and 
measurement  noise.  This  analysis  is  consistent  with  related  research  which  showed  that 
shot  noise  is  correlated  with  respect  to  spatial  frequency  [54].  The  amount  of  correlation 
depends  on  the  product  of  the  mean  OTF  and  the  mean  object  spectrum  at  a  difference 
frequency.  This  linear  filter  can  provide  a  useful  alternative  to  nonlinear  iterative  techniques 
when  appropriate  statistical  models  are  available. 

Chapter  IV  is  organized  as  follows.  Section  4.2  outlines  the  complete  derivation  of 
the  new  filter  for  both  random  and  deterministic  OTF  cases.  Section  4.3  gives  a  signal-to- 
noise  ratio  (SNR)  interpretation  to  the  estimator.  In  Section  4.4,  theoretical  filter  mean 
square  error  (MSE)  expressions  are  derived  for  the  vector  Wiener  filter,  scalar  Wiener  filter, 
and  the  unfiltered  data.  These  expressions  are  important  for  demonstration  and  analysis  in 
Chapters  V  and  VI.  Section  4.5  provides  alternate  filter  expressions  which  incorporate  the 
mean  and  covariance  of  the  random  quantities.  Finally,  some  comments  are  made  in  Section 
4.6  regarding  the  optimality  of  the  new  filter  expressions.  The  chapter  ends  with  a  brief 
summary  in  Section  4.7. 

4.2  Fourier  Domain  Filter  Derivation 

In  Chapter  II,  the  potential  performance  advantages  associated  with  the  vector  Wiener 
filter  were  noted  for  non-stationary  image  ensembles.  However,  short  exposure  images  col¬ 
lected  through  atmospheric  turbulence  are  blurred  by  an  unknown  random  OTF.  In  addition, 
measurement  noise  further  degrades  optical  system  resolution.  Equation  (2.36)  does  not  ac¬ 
count  for  these  factors  and  must  be  extended  for  this  imaging  application. 
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4.2.1  Wiener-Hopf  Expression.  Let  us  consider  an  N  x  N  image  array.  Based  on 
the  Fourier  domain  vector-matrix  image  model  expression  given  in  Eq.  (2.9),  the  objective 
is  to  find  an  N'^  x  Fourier  domain  filter  matrix  Mr  such  that 


6  =  MrD, 


(4.1) 


where  the  MSE 

e*  =  E  [(O  -  6)"(o  -  6)] 


(4.2) 


is  minimized.  Using  the  matrix  trace  operator  Tr{*}  [35],  Eq.  (4.2)  can  be  written  as 


=  E  [Tt  {(O  -  6)(0  -  6)"}]  . 


(4.3) 


Substituting  Eq.  (4.1)  into  Eq.  (4.3),  bringing  the  expectation  operator  inside  the  trace 
operator,  and  expanding  the  expression  yields 


=  Tr  |i?oo  -  MrRqd  ~  RodMr  +  MrRddMr  | ,  (4.4) 


where  Rqo  =  E[00^]  is  the  object  Fourier  domain  autocorrelation  matrix.  Rod  =  E[OD^] 
is  the  object-detected  image  Fourier  domain  cross-correlation  matrix,  and  Rdd  =  E[DD^]  is 
the  detected  image  Fourier  domain  autocorrelation  matrix.  To  find  the  filter  transformation 
matrix  Mr  that  minimizes  take  the  derivative  of  Eq.  (4.4)  with  respect  to  Mr,  set  this 
derivative  equal  to  the  zero  matrix,  and  solve  for  Mr.  The  resultant  derivative  is  [35] 


de^ 

BMr 


—  —2Rod  +  2i?£)£)M^  =  0. 


(4.6) 


Thus,  the  linear  minimum  MSE  Fourier  spectrum  estimate  is 


6  —  Rod{Rdd) 


(4.6) 


where  the  transformation  Mr  =  Rod{Rdd)~^  satisfies  a  Wiener-Hopf  equation  [43]. 


4-2 


4.2.2  Object- Detected  Image  Cross-correlation.  To  write  the  new  vector  Wiener 
filter  expression,  the  Fourier  domain  correlation  matrices  Rqd  and  Rdd  must  be  derived 
using  the  image  model  given  in  Eq.  (2.7).  First,  consider  the  cross-correlation  between  the 
object  spatial  frequency  (u,  v)  and  detected  image  spatial  frequency  («',  v')  denoted 


Roniu,  v;  u\  v')  =  E[0{tf,  v)D*{u',  ?;')]. 


(4.7) 


Substituting  Eq.  (2.7)  into  Eq.  (4.7)  and  writing  as  the  sum  of  two  terms  yields 


Rod{u,v,u',v')  =  E 


K 

I 

n=l 


Oiu,  v)  exp  {j2TT{u'Xn  +  v'yn)} 


+  E 


0(u,  v)  Y  %  exp  {j2Tr{u'xp  +  v'yp)} 

p=i 


(4.8) 


The  second  term  in  Eq.  (4.8)  is  zero  because  Up  is  both  independent  of  O  and  zero  mean. 
The  first  term  can  be  evaluated  using  nested  conditional  expectations  following  the  technique 
presented  in  Refs.  [23]  and  [74].  The  random  quantities  are  the  object  spectrum,  the  OTF, 
the  total  number  of  photoevents  K,  and  the  photon  arrival  location  {xn,yn)‘  Bayes  rule  [62] 
can  be  used  to  rewrite  the  joint  probability  density  function  (PDF)  in  the  first  term  using 
conditional  PDFs.  These  conditional  PDFs  translate  to  conditional  expectations  such  that 


Rod{u,V]u',v')  =  Eo[E7.^|o(Eif|7^,o{Ea;„,y„|ir,W,o[0(M,  f) 

K 

X  Y  6xp  {j27r(u  'x„  +  t;'yn)}]})]  ,  (4.9) 

n=l 

where  the  notation  Ea\b  denotes  the  conditional  expected  value  of  the  random  event  A 
given  B.  An  expression  similar  to  the  innermost  conditional  expectation  of  Eq.  (4.9)  has 
been  evaluated  previously  [74].  The  derivation  details  in  this  case  are  given  in  Appendix  A.l 
with  the  final  result  written  as 


E, 


Xn,yn\K,'H,0 


K 


v)  Y  exp  {j2TT{u'Xn  +  v'yn)} 

n~l 


v')0{u,  v)0*n{u',  v'),  (4.10) 
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where  the  normalized  object  Fourier  spectrum  0„(«,  v)  is  defined  as 


On{u,v) 


0{u,v) 

0(0,0) 


0{u,v) 

K 


(4.11) 


Evaluation  of  the  remaining  nested  expectations  is  trivial  since  K,  H,  and  O  are  all  mutually 
independent.  This  mutual  independence  converts  all  the  nested  quantities  into  an  uncoupled 
product  of  expectations,  which  gives  the  following  expression  for  Eq.  (4.7) 


Rod{u,  v]  u',  v')  =  K  H\u',  u')E  [0(m,  v)Ol{u',  ?;')] , 


(4.12) 


where  'H(u,  v)  is  the  mean  OTF.  Clearly,  the  expectation  above  is  the  object  Fourier  domain 
autocorrelation,  except  for  a  normalization  factor  associated  with  0(u,v).  Thus,  Eq.  (4.11) 
can  be  used  to  write  Eq.  (4.12)  as 


Kn’(u',v')ElO(u,v)0:(u',v')j  =  (K)^  ■H'(n’,v')E[0„(u,v)0:(nW)] 


(4.13) 


Finally,  Eq.  (4.7)  becomes 

Eod(u,v;u',v')  =  (Kf  7{*(u',v')RonO„(u,v;u',v'),  (4.14) 

where  RonOni'^^ autocorrelation  between  the  (u,  v)  and  the  {u',  v')  spatial  fre¬ 
quencies  of  the  normalized  object  Fourier  spectrum.  Equation  (4.14)  can  be  expressed  using 
the  vector-matrix  notation  introduced  in  Chapter  II.  Note  that  the  functional  dependence 
of  the  OTF  on  the  spatial  frequency  (u',  v')  is  equivalent  to  multiplying  a  particular  column 
of  the  matrix  i2o„o„  by  H*  {u' ,  v').  This  process  is  equivalent  to  multiplying  by  the  diagonal 
matrix  such  that  the  final  result  for  the  object-detected  image  Fourier  domain  correlation 
matrix  is 

Rod  =  {KY  RonOn'Rd-  (4-15) 
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4.2.3  Detected  Image  Autocorrelation.  To  complete  the  derivation,  the  correlation 
between  detected  image  spatial  frequencies  (w,  v)  and  (u',  v')  denoted 


Rdd{u,  v\ u',  v')  =  E  [D{u,  v)D*{u',  v')] , 


(4.16) 


must  be  derived.  Substituting  the  image  model  given  in  Eq.  (2.7)  into  Eq.  (4.16)  and 
expanding  the  expression  within  the  expectation  yields 


Rdd{u,v]u',v')  =  E 


■  K  K 

Y.  Y  {-j27r(MrCn  -  u'xm  +  vyn  -  v'ym)} 

.71=1  m=l 
p  p 

exp  {-j27r(uXp  -  u'xq  +  vyp  -  v'yq)} 

p=l  g=l 

Rdd^u,  V, »')  +  Rooi'^, »'). 


(4,17) 


where  the  cross  terms  are  zero  because  Up  is  zero  mean  and  independent  of  K,  and  O. 
Consider  the  first  term,  and  write  using  nested  conditional  expectations 

[23, 74]  as  before  which  yields 


=  ^o\^'H\o{^K\H,o\^Xn,yn,Xm.ym\K,H,0 

K  K 

X  ^  exp  {-j2'K{uXn  -  u'xm  +  vyn  “  v'ym)}] })] .  (4.18) 

71=1  m=l 

The  double  summation  in  ^xn,yn,xm,ym\K,'H,o  has  two  types  of  terms:  K  terms  in  which  n  =  m 
and  -  K  terms  in  which  n  ^  m.  Thus,  'Ea;„,y„,xm,ym\K,n,o  in  Eq.  (4.18)  can  be  rewritten 
in  a  different  form: 


E, 


Xn  ^yn  jXm  jym 


r  K 


K  K 


Y  Y  {-j^T^iuXn  -  u'Xm  +  Vyn  “  v'ym)} 

n=l  m=l 


,n=l 


^Xn,yn\K,'H,0 

+  ^Xn,yn,Xm,ym\K,'H,0 


Y  ®xp  {-j27r((u  -  u')Xn  +  {v-  v')ym.)} 


K  K 


Y  Y  ®XP  {-j27r(ua;„  -  u'xm  +  vy^  -  v'ym)} 

,71=1  771=1 


.(4.19) 


n^m 
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The  derivation  details  associated  with  Eq.  (4.19)  are  given  in  Appendix  A. 2  with  the  final 
result  written  as 


W  =  (K^-K)0„(u,v)0:(u',v')n(u,v)n'(u',v') 

+  K'H{u  —  u',v  —  v')On{u  —  u' ,v  —  v'),  (4.20) 

where  the  •  notation  represents  the  bracketed  quantity  on  the  left  side  of  Eq.  (4.19).  With 
Eq.  (4.20)  in  hand,  the  remaining  conditional  expectations  can  be  evaluated  as  before  by 
noting  the  mutual  independence  of  K,  H,  and  O,  which  leads  to 

■Rdd(«>  RonOni'f^,  u',  v')Rhh{u,  V,  u',  v') 

+  K  7{{u  —  u',v  —  v')On(u  —  u',v  —  v'),  (4-21) 

where  V]  u',  v')  is  the  autocorrelation  between  (u,  v)  and  [u' ,  v')  spatial  frequencies  of 

the  OTF  and  On{u^v)  denotes  the  normalized  mean  object  spectrum.  The  random  variable 
K  is  conditionally  Poisson  distributed,  given  H  and  O  [74].  Therefore,  the  second  moment 
of  K  can  be  written  as  [62, 74] 


^K‘^]  =K^  =  K  +  {Kf.  (4.22) 

Using  Eq.  (4.22),  Eq.  (4.21)  can  be  written  in  its  final  form: 

=  {KyRonOn(u,v;u',v')Rnn{u,v,u',v') 

+  K  H{u  —  u',v  —  v')On{u  —  u',v  —  v').  (4.23) 

The  second  term  in  Eq.  (4.17)  was  evaluated  previously  [74]  for  the  assumptions  given  in 
Chapter  II  such  that 

w',  v')  =  Pul5{u  -u',v  -  v'),  (4.24) 
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where  P  =  is  the  number  of  pixels  in  the  detector.  Combining  the  results  from  Eqs.  (4.23) 

and  (4.24)  gives  the  final  expression  for  R£,d{u,v;  u',v'): 

Rdd{u,v]u',v')  =  (K)^Ro^On{u,v;u',v')Rn'H{u,v,u',v') 

+  K  H{u  —  u',V  —  v')On{u  —  u',v  —  v') 

+  Pol5{u  —  u\v  —  v').  (4.25) 

Equation  (4.25)  can  be  written  conveniently  using  vector-matrix  notation  as  before,  which 
yields  a  compact  expression  for  the  matrix  Rod 

Rod  =  {^Y^OnOn  ®  Rhh  +  Rnni  (4.26) 

where  the  measurement  noise  Fourier  domain  autocorrelation  between  (u,  u)  and  {u',v') 
spatial  frequencies  is  defined  as 


Rnn{u,  v]  u',  v')  =  K  H{u  -  u',v  -  v')On{u  -  u',v  -  v')  +  Pal5{u  -  u',v  -  v').  (4.27) 

4.2.4  Final  Result.  Returning  to  Eq.  (4.6),  inserting  the  results  given  in  Eqs.  (4.15) 
and  (4.26),  and  dividing  through  by  the  scalar  {K)"^  give  the  final  expression  for  the  extended 
vector  Wiener  filter  [14]: 


O  —  RonOnFtd  {RouOn  ®  Rhh  +  K  Rnn]  D>  (4.28) 

where  the  elements  of  the  matrix  iZivjv  are  defined  in  Eq.  (4.27).  Equation  (4.28)  is  the 
main  result  of  this  section  and  the  subject  of  discussion  and  experimentation  in  Chapters 
V  and  VI.  Equation  (4.28)  is  new  and  has  only  recently  been  applied  to  images  degraded 
by  a  random  OTF  and  measurement  noise  [12-14].  If  the  OTF  is  deterministic  and  known, 
Eq.  (4.28)  can  be  written  in  the  following  form: 

6  =  iJo.o.W;  {HiRo.o,K  +  D,  (4.29) 
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where  the  deterministic  OTF  is  substituted  for  the  mean  OTF  when  defining  the  Fourier 
domain  noise  autocorrelation  matrix  Rnn- 


4.3  Signal-to-Noise  Ratio  Interpretation 

It  is  convenient  to  express  the  measurement  noise  Fourier  domain  correlation  in  terms 
of  SNR  quantities.  Multiplying  Rnn  by  the  scalar  K  as  required  by  Eq.  (4.28)  yields 

K  ‘^Ri^r^f(u,v,u',v')  =  K  ^H{u  — u',v  —  v')On{u  —  u',v  —  v') 

+  ^Pa^/{K)^'^6{u  -  u',v  -  v').  (4.30) 


The  single  frame  image  spectrum  SNR  is  defined  as  [74] 


SNRi(w,u)  = 


E[l£>(«,^;) 


{Var[D(tf,  t))]}^/2’ 


(4.31) 


and  was  calculated  previously  for  an  image  degraded  by  measurement  noise,  yielding  [74] 


SNRi(m,u)  = 


\H{u,  v)\\On{u,v)\ 


{\On(u,v)\^Y8.T[niu,v)]  +  K~' +Pa^J{Kyy/^‘ 


(4.32) 


where  Var[«]  denotes  the  variance  of  the  bracketed  expression.  The  last  two  quantities  in  the 
denominator  of  Eq.  (4.32)  are  related  to  the  shot  noise  and  the  signal-independent  detector 
read  noise,  respectively.  Consider  the  following  SNR  expressions  [90] 


SNRfc  =  ^/k,  (4.33) 

SNR,  =  -^=,  (4.34) 

where  SNRfc  is  associated  with  the  shot  noise  and  SNR,  is  associated  with  the  detector  read 
noise.  Equations  (4.33)  and  (4.34)  were  used  previously  to  rewrite  the  single  frame  image 
spectrum  SNR  [90] .  In  the  same  way,  the  measurement  noise  Fourier  domain  correlation  in 
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Eq.  (4.28)  can  be  written  as 


K  '^Rnn{u,v,u',v')  =  SNR^^H{u  —  u',v  —  v')Oniu  —  u' ,v  —  v') 

+  S^R;^5{u-u',v-v').  (4.35) 

These  SNR  quantities  control  the  amount  of  regularization  provided  to  each  spatial  frequency 
correlation  component  by  the  vector  Wiener  filter.  When  the  measurement  noise  is  low  (i.e., 
high  light  level  and  minimal  detector  read  noise),  negligible  regularization  occurs  and  the 
filter  is  faithful  to  the  data.  When  the  measurement  noise  is  high  (i.e.,  low  light  level 
or  significant  detector  read  noise),  regularization  is  applied  based  on  the  spatial  frequency 
correlation  of  the  shot  noise  and  the  strength  of  the  uncorrelated  detector  read  noise.  Under 
extremely  high  noise  conditions,  the  vector  Wiener  filter  will  smooth  the  data,  but  always 
consistent  with  the  known  object  Fourier  domain  statistical  model.  In  the  next  section, 
MSE  expressions  are  introduced  which  will  be  important  for  filter  performance  comparison 
in  Chapters  V  and  VI. 

4.4  Mean  Square  Error  Expressions 

The  vector  Wiener  filter  given  in  Eq.  (4.28)  minimizes  the  scalar  MSE  as  given  in 
Eq.  (4.2).  Thus,  filter  performance  can  be  analyzed  by  examining  the  statistics  of  the  zero 
mean  error  between  the  true  and  estimated  Fourier  spectra  e  =  O  —  6.  The  error  correlation 
matrix  associated  with  a  Fourier  domain  linear  minimum  MSE  filter  can  be  written  as  [43] 

R,,  =  E  [(O  -  6)(0  -  6)^]  .  (4.36) 

In  Eq.  (4.1),  a  linear  estimator  model  was  considered  which  incorporated  a  specific  vector 
Wiener  filter  transformation  matrix  Mr.  Here,  the  same  linear  model  will  be  used  but  with 
an  arbitrary  transformation  matrix  Mx  such  that  Eq.  (4.1)  becomes 


6  =  MxD. 


(4.37) 
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Substituting  Eq.  (4.37)  for  6,  expanding  terms,  and  evaluating  the  expectation  yields  [14,43] 
Ret  =  Roo  —  ^xR-dRoO  —  RooR-d^X  +  ^X  {RoO  ©  RnH  +  Rnn)  Mx-  (4.38) 
For  the  vector  Wiener  filter  transformation  Mr  given  in  Eq.  (4.28),  Eq.  (4.38)  reduces  to 
iiee  =  Roo  —  MpTidRoo  (Vector  Wiener  Filter).  (4.39) 


A  similar  expression  is  also  available  for  the  analogous  scalar  Wiener  filter  whose  diagonal 
transformation  matrix  Ms  has  non-zero  elements  defined  as 


Ms{u,v) 


(u,  v) 

u)|2]  +  (k  +  Pcr2)  /g,(u,  v)  ’ 


(4.40) 


where  Go{u,v)  is  the  PSD  of  the  object  as  used  in  Eq.  (2.35).  The  error  correlation  matrix 
associated  with  the  scalar  Wiener  filter  is 


Ree  =  Roo  —  MsRdRoo  ~  RooRd^s  +  {Roo  ©  Rhh  +  Rnn) 

(Scalar  Wiener  Filter).  (4.41) 

Finally,  it  is  often  of  interest  to  know  the  baseline  error  statistics  associated  with  the  un¬ 
filtered  detected  image.  When  Mx  =  I,  where  I  is  the  identity  matrix,  Eq.  (4.37)  gives 
6  =  D  as  required  and  the  error  correlation  matrix  becomes 

Ree  =  Roo  —  RdRoo  ~  RooRd  +  {Roo  ©  Rhh  +  Rnn) 

(Detected  Image).  (4.42) 

Equations  (4.39),  (4.41),  and  (4.42)  are  only  dependent  on  the  statistics  of  the  object  class, 
OTF,  and  noise.  Thus,  a  theoretical  performance  study  based  on  these  error  statistics  is 
possible  without  running  a  Monte  Carlo  simulation.  In  the  next  section,  two  alternate  vector 
Wiener  filter  expressions  are  derived  which  incorporate  mean  and  covariance  statistics. 
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4-5  Alternate  Filter  Expressions 

The  vector  Wiener  filter  given  in  Eqs.  (4.28)  and  (4.29)  is  a  function  of  correlation 
matrices.  In  many  statistical  signal  processing  and  image  processing  texts,  Bayesian  linear 
minimum  MSE  filters  are  expressed  in  terms  of  mean  and  covariance  [41,43].  Thus,  two 
complementary  vector  Wiener  filter  equations  are  given  below  in  the  interest  of  completeness. 

Let  us  now  reconsider  the  form  of  the  random  object  realization  O.  Such  a  random 
vector  can  always  be  written  as  a  function  of  its  mean  O  and  a  random,  zero  mean  component 
AO  yielding 

0  =  0  +  AO,  (4.43) 

where  AO  is  zero  mean  with  covariance  matrix  Cqo-  In  general,  a  complex  covariance 
matrix  is  defined  as  [62] 

Cxx  =  E  [(X  -  X)(X  -  X)"^]  ,  (4.44) 

where  X  is  a  complex  random  vector  with  mean  X.  The  object  estimate  vector  can  now  be 
written  as 


6  =  o  +  a6 

=  O  +  McAD,  (4.45) 

where  AD  =  D  —  D,  D  is  the  mean  detected  image  vector,  and  Me  is  a  new  filter  trans¬ 
formation  matrix.  Following  analysis  similar  to  that  given  in  Eqs.  (4.1)-(4.5)  produces  a 
Wiener-Hopf  expression  incorporating  covariance  matrices 


AO  —  Cod  {Cdd)  ^  (4.46) 

where  Cqd  is  the  cross-covariance  between  object  and  detected  image  and  Cdd  is  the  au¬ 
tocovariance  of  the  detected  image.  The  analysis  in  Section  4.2  will  now  be  used  to  find 
expressions  for  Cqd  and  Cdd- 
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By  definition,  the  matrix  Cqd  can  be  written  as  [62] 


Cod  =  Rod-0  (4.47) 

The  matrix  Rod  was  derived  above  and  is  given  in  Eq.  (4.15).  Previous  analysis  has  shown 
that  the  mean  detected  image  vector  is  D  =  iif  'HdQn  [74].  Substituting  Eq.  (4.15)  and  the 
mean  detected  image  Fourier  spectrum  into  Eq.  (4.47)  yields 

Cod  =  (Kflto„oM-(Kyo„o”% 

=  (jSo.o.  -  0„  o")  TT, 

=  (KfCo^oM  (4.48) 

The  matrix  Cdd  can  also  be  written  in  a  form  which  incorporates  Rod  and  D  such  that  [62] 

Cod  —  Rdd  —  D  D  .  (4.49) 

Rod  was  derived  above  and  is  given  in  Eq.  (4.26).  Substituting  Eq.  (4.26)  and  the  mean 
detected  image  Fourier  spectrum  into  Eq.  (4.49)  yields 

Cdd  =  {KyRonOn  ®  ^hh  +  Rnn  —  {RY  (o„  ©  (h  ,  (4.50) 

where  H  denotes  a  P-length  vector  containing  the  mean  OTF  elements  as  noted  in  Chapter 
II.  Equation  (4.50)  can  be  rewritten  in  terms  of  covariance  matrices  using  the  general  form 
shown  in  Eqs.  (4.47)  and  (4.49).  Combining  Eq.  (4.48)  with  the  modified  form  of  Eq.  (4.50) 
gives  an  alternate  filter  expression  incorporating  covariance  matrices 

6  =  KOn^Co.ojl*i 

X  {ConOn  0  Cun  +  Cnn  0  (On  +  C7o„o„  0  +  R  “^Cnn) 

X  (p-RTtdOn),  (4.51) 
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where  Cnn  =  Rnn  since  the  measurement  noise  is  zero  mean  and  is  defined  as  given 
in  Eq.  (4.27).  When  the  OTF  is  deterministic  and  known,  Eq.  (4.51)  becomes 

6  =  KOn  +  Co^oM  {WonoM  +  (d  -  ¥  HdOn)  .  (4.52) 

The  next  section  provides  some  brief  comments  related  to  the  optimality  of  the  new  vector 
Wiener  filter.  Here,  the  linear  constraint  will  be  investigated  with  respect  to  the  underlying 
PDFs  associated  with  object  and  noise. 

4.6  Comments  on  Filter  Optimality 

The  vector  Wiener  filter  is  optimal  with  respect  to  MSE  or  error  variance,  regardless  of 
the  stationarity  of  the  underlying  random  processes.  The  scalar  Wiener  filter  is  only  optimal 
in  this  sense  when  object  and  noise  are  stationary.  Here,  an  “optimal”  filter  is  one  that 
provides  the  minimum  error  variance  for  a  particular  imaging  scenario  and  filter  class.  The 
filter  class  referred  to  here  is  the  class  of  Bayesian  linear  filters  only.  A  nonlinear  estimator 
may  exist  that  provides  an  ensemble  of  solutions  with  less  error  variance  than  the  vector 
Wiener  filter  as  illustrated  in  Fig.  4.1.  In  the  context  of  this  figure,  the  vector  Wiener  filter 
is  guaranteed  to  provide  minimum  error  variance  among  the  class  of  linear  Bayesian  filters 
represented  by  the  small  circle  marked  “Linear” ,  not  the  larger  circle  marked  “General” . 

As  noted  in  Chapter  II,  linear  minimum  MSE  estimators  have  been  extensively  studied 
and  applied  to  a  wide  variety  of  signal  processing  problems  [43].  This  body  of  knowledge 
reveals  one  important  case  in  which  the  linear  estimator  is  optimal  with  respect  to  MSE  for 
the  class  of  all  Bayesian  filters.  This  situation  exists  when  the  noise  and  object  are  Gaussian 
distributed  and  statistically  independent  [43].  However,  this  case  does  not  apply  to  the 
image  model  given  in  Eq.  (2.6)  since  the  shot  noise  is  not  independent  of  the  signal.  The 
rest  of  this  section  will  investigate  the  effect  of  the  semi-classical  model  on  the  optimality 
of  the  vector  Wiener  filter.  Section  4.6.1  will  review  key  background  details  related  to  the 
linear  Bayesian  estimator.  Section  4.6.2  presents  the  signal-independent  noise  case.  Finally, 
Section  4.6.3  outlines  optimality  issues  related  to  the  semi-classical  image  model. 
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Vector  Wiener  Filter 


Figure  4.1  The  vector  Wiener  filter  is  optimal  with  respect  to  minimum  error  variance 
only  among  the  class  of  linear  filters  represented  by  the  small  circle  marked 
“Linear”.  Further  statements  with  respect  to  optimality  require  assumptions 
regarding  object  and  data  distributions. 
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4.6.1  Background.  By  definition,  the  Bayesian  minimum  MSB  estimator  is  the 
mean  of  the  posterior  PDF,  represented  by  Eo|d[0|D]  [43].  If  Eo|d[0|D]  can  be  shown  to 
be  linearly  related  to  the  data,  then  the  vector  Wiener  filter  must  be  the  minimum  error 
variance  filter  over  the  class  of  all  Fourier  domain  Bayesian  minimum  MSE  filters.  This 
special  case  is  analogous  to  constraining  the  “best”  estimator  to  exist  within  the  small 
circle  marked  “Linear”  in  Fig.  4.1.  Thus,  to  understand  these  optimality  issues  fully,  the 
conditions  under  which  Eo|d[0|D]  is  linearly  related  to  the  data  must  be  investigated.  To 
lay  the  groundwork  for  this  investigation,  consider  two  jointly  complex  Gaussian  random 
vectors  X  and  Y  whose  vector  [X^Y^]^  can  be  described  as 

X  r  ^  ^XX  CxY 

~  AT  _  ,  ^  (4  53) 

Y  \  ^  ^Yx  Cyy  j 

where  the  notation  M  denotes  a  random  quantity  with  Gaussian  distribution,  given  mean 
vector,  and  given  covariance  matrix.  An  equivalent  representation  of  Eq.  (4.53)  is  that  the 
joint  PDF  px,y(X,  Y)  is  Gaussian.  Since  a  Gaussian  PDF,  real  or  complex,  is  completely 
described  by  its  mean  and  covariance,  Eq.  (4.53)  provides  a  complete  statistical  description. 
Further,  it  has  been  shown  that  if  px,y(X,  Y)  is  complex  Gaussian  as  described  in  Eq.  (4.53), 
the  conditional  PDF  py|x(Y|X)  is  also  complex  Gaussian  with  mean  [43] 

Ey|x[Y|X]  =  Y  +  Cyx  {CxxT"  (X  -  X).  (4.54) 

Equation  (4.54)  is  the  mean  of  the  posterior  PDF  and  also  the  Bayesian  minimum  MSE 
estimator  of  Y.  Thus,  if  X  and  Y  are  jointly  complex  Gaussian,  Ey|x[Y|X]  is  linearly 
related  to  X.  When  the  vector  D  is  substituted  for  X  and  the  vector  O  is  substituted  for 
Y,  Eq.  (4.54)  can  be  viewed  in  terms  of  the  image  model  given  in  Eq.  (2.9),  repeated  as 


D  =  W  0  O  +  N, 


(4.55) 


and  the  vector  Wiener  filter  given  in  Eq.  (4.52).  Therefore,  if  D  and  O  are  jointly  complex 
Gaussian,  the  general  Bayesian  minimum  MSE  estimate  is  linear  to  the  data  and  the  vector 
Wiener  filter  is  the  general  minimum  error  variance  estimator. 

Under  what  conditions  are  D  and  O  jointly  complex  Gaussian?  The  following  dis¬ 
cussion  will  answer  this  question  based  on  statistical  assumptions  about  the  general  linear 
model  given  in  Eq.  (4.55). 

4.6.2  Signal-Independent  Noise.  Let  us  take  a  step  back  from  the  image  degrada¬ 
tion  model  first  given  in  Chapter  II  and  assume  that  the  OTF  is  deterministic.  In  addition, 
let  us  assume  that  the  complex  noise  vector  N  is  independent  of  O  and  Gaussian  distributed. 
This  assumption  is  used  in  many  image  processing  applications  [41]  but  is  not  consistent  with 
an  image  degraded  by  shot  noise  or  the  image  model  in  Eq.  (2.6).  The  joint  PDF  ,o(D,0) 
can  be  rewritten  in  terms  of  conditional  and  marginal  PDFs  such  that  [62] 

Pd,o(D)0)  =f>D|o(D|0)po(0)-  (4.56) 

Since  the  noise  has  a  Gaussian  PDF  and  fC  is  deterministic,  Pd|o(D]0)  is  also  Gaussian 
distributed.  Thus,  if  po(0)  is  Gaussian  distributed,  the  joint  PDF  is  Gaussian  [62]  and  the 
vector  Wiener  filter  is  the  general  Bayesian  minimum  error  variance  filter. 

4.6.3  Semi- Classical  Model.  While  the  Gaussian  distributed  signal-independent 
noise  case  is  valid  in  many  applications,  it  does  not  adequately  represent  low  light  imaging 
of  atmospheric  turbulence-degraded  images.  Here,  shot  noise  and  the  randomness  of  the 
OTF  introduce  dependence  on  the  signal.  To  consider  this  case,  it  is  helpful  to  rewrite  the 
OTF  in  the  following  form 

Hd  =  nd  +  AUd,  (4.57) 

where  AHd  is  a  diagonal  matrix  representing  a  realization  of  the  random  OTF  component. 
Now  Eq.  (4.55)  can  be  rewritten  using  Eq.  (4.57)  such  that 

D  =  HdO  -f  Nr,  (4.58) 
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where  the  total  noise  is  defined  as  Ny  =  AHdO  +  N  and  is  dependent  on  the  object  via 
the  random  OTF  and  shot  noise  components.  Even  if  the  object  PDF  po{0)  is  Gaussian 
distributed,  the  joint  PDF  Pd,o(D,  O)  will  not  be  Gaussian  since  Ny  is  not  strictly  Gaussian 
and  is  dependent  on  the  object.  Further,  the  signal  dependence  of  the  noise  precludes  the 
existence  of  any  case  where  Pd,o(D,0)  is  Gaussian  distributed  regardless  of  the  choice  of 
noise  or  object  PDF.  Thus,  in  a  strict  sense,  no  statements  about  the  optimality  of  the  vector 
Wiener  filter  can  be  made  outside  of  the  class  of  Bayesian  linear  filters. 

However,  in  a  limiting  sense,  the  vector  Wiener  filter  approaches  optimality  across  the 
class  of  all  Bayesian  filters  when  two  conditions  are  present:  the  average  photon  count  K 
grows  large  and  the  OTF  covariance  grows  small.  In  Appendix  B,  it  is  shown  that  pn  |o(D|0) 
approaches  a  complex  Gaussian  distribution  when  the  total  photon  count  K  is  sufficiently 
large,  via  a  Central  Limit  Theorem  argument  [23,62].  For  most  practical  images,  K  is  very 
large,  on  the  order  of  1,000  to  1,000,000  photoevents.  This  value  is  more  than  adequate  for 
valid  application  of  the  Central  Limit  Theorem  [62].  The  additional  assumption  of  small 
OTF  covariance  supports  the  argument  that  Nr  approaches  a  Gaussian  distribution.  In 
addition,  these  two  conditions  yield  the  signal-independent  noise  case.  Thus,  when  po(0) 
is  Gaussian,  the  joint  PDF  approaches  a  Gaussian  distribution,  and  the  vector  Wiener  filter 
can  be  viewed  as  essentially  optimal  for  all  Bayesian  minimum  error  variance  filters.  The 
two  conditions  noted  above  represent  imaging  under  high  light  level  conditions  where  an 
accurate  estimate  of  the  turbulence-induced  OTF  is  available. 

4.7  Summary 

In  this  chapter,  a  Fourier  domain  filter  was  derived  which  incorporates  object,  blur,  and 
noise  statistical  models.  This  analysis  extends  the  vector  Wiener  filter  to  account  properly  for 
both  shot  noise  and  detector  read  noise  as  modeled  in  Eq.  (2.7).  The  shot  noise  correlation 
depends  on  the  product  of  the  mean  OTF  and  the  mean  object  spectrum  at  a  difference 
frequency  [14,54].  The  theoretical  optimality  of  the  vector  Wiener  filter  was  also  investigated 
with  respect  to  assumptions  about  the  underlying  distributions  of  the  data  and  object.  When 
the  detected  image  is  degraded  by  signal-dependent  shot  noise  and  a  random  OTF,  the  vector 
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Wiener  filter  is  only  guaranteed  to  be  the  minimum  error  variance  estimator  among  the  class 
of  linear  Bayesian  filters. 

In  the  next  chapter,  simulated  binary  star  data  is  processed  with  the  vector  Wiener 
filter.  The  binary  star  represents  an  important  astronomical  imaging  application.  These 
results  include  data  degraded  by  both  a  fixed  and  a  random  OTF. 


4-18 


V.  Vector  Wiener  Filter  Processing  of  Binary  Star  Pairs 

5.1  Introduction 

In  the  previous  chapter,  a  Fourier  domain  vector  Wiener  filter  was  derived  which  in¬ 
corporates  complete  object,  blur,  and  noise  correlation  statistics.  Equation  (2.7)  was  used  in 
the  derivation  to  model  all  degradation  effects  in  the  detected  image.  This  model  is  consis¬ 
tent  with  an  image  degraded  by  atmospheric  turbulence,  shot  noise,  and  detector  read  noise. 
In  this  chapter,  the  vector  Wiener  filter  is  used  to  process  an  important  non-stationary 
astronomical  object  class,  the  binary  star  pair.  The  goal  is  to  study  vector  Wiener  filter 
performance  when  the  object  and  optical  transfer  function  (OTF)  statistics  are  known  ex¬ 
actly.  Here,  performance  is  compared  against  the  scalar  Wiener  filter  using  mean  square 
error  (MSE),  a  correlation  coefficient,  mean  square  phase  error  (MSPE),  and  image  real¬ 
izations  for  both  a  deterministic  and  random  OTF.  In  all  cases,  the  vector  Wiener  filter 
provides  reconstructions  that  are  superior  to  those  of  the  scalar  Wiener  filter.  The  results 
also  illustrate  the  superresolution  capability  of  the  vector  Wiener  filter.  In  this  dissertation, 
superresolution  is  defined  as  the  extension  of  the  detected  image  Fourier  spectrum  to  regions 
where  no  data  was  measured  [55] . 

Chapter  V  is  organized  as  follows.  Section  5.2  provides  general  details  of  the  Monte 
Carlo  simulation  used  to  generate  the  data  to  include  information  about  the  random  object, 
measurement  noise,  and  performance  metrics.  Sections  5.3  and  5.4  present  the  data  and 
conclusions  associated  with  the  deterministic  and  random  OTF  cases,  respectively.  The 
chapter  ends  with  a  summary  in  Section  5.5. 

5.2  Simulation 

In  this  section,  details  of  a  simulation  used  to  study  vector  Wiener  filter  performance 
are  provided.  The  term  simulation  is  used  to  refer  to  Monte  Carlo  experiments  involving  the 
repeated  application  of  the  filter  to  random  draws  of  object,  OTF,  and  measurement  noise. 
The  discussion  here  includes  details  associated  with  random  object  and  measurement  noise. 
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5.2.1  Random  Object.  To  compare  vector  and  scalar  Wiener  filter  performance, 
a  class  of  simple  non-stationary  random  objects  is  needed  with  Fourier  domain  statistics 
that  can  be  derived  analytically.  The  data  generated  for  this  study  are  based  on  a  pair  of 
Gaussian  functions  with  locations  that  are  random.  Figure  5.1  shows  (a)  16  x  16  and  (b) 
32  X  32  detector  array  object  realizations  that  can  be  expressed  mathematically  as  [22] 


o{x,  y)  =  hp  exp  <  — tt 


'{x-XpY +  {y-ypy 


Wr, 


+  hs  exp  <  — TT 


'(x  -  Xsf  +  {y-  ysf 


w. 


(5.1) 


where  {Xp,  yp)  is  the  location  of  the  primary  function,  {xg,  yg)  is  the  location  of  the  secondary 
function,  hp  is  the  peak  irradiance  of  the  primary,  hg  is  the  peak  irradiance  of  the  secondary, 
Wp  is  the  primary  width  parameter,  and  Wg  is  the  secondary  width  parameter.  The  func¬ 
tion  locations  (xp,yp)  and  {xs,ys)  are  independent,  uniformly  distributed  random  vectors 
restricted  to  a  W  by  W  pixel  region  in  the  center  of  the  image  plane.  For  example,  W  =  1 
requires  both  functions  to  be  located  at  the  center  of  the  image  plane,  and  the  object  realiza¬ 
tion  is  deterministic.  Since  the  function  locations  have  a  uniform  distribution,  the  parameter 
W  can  also  be  viewed  as  the  dimension  of  a  spatial  domain  support  constraint.  As  shown 
in  Fig.  5.1,  the  normalized  peak  irradiance  of  the  primary  is  hj,  =  1  and  the  secondary  is 
hg  =  0.5.  Both  primary  and  secondary  functions  have  width  parameters  Wp  =  Wg  =  0.5 
pixels  for  the  16  x  16  detector  array  in  (a)  and  Wp  —  Wg  =  1  pixel  for  the  32  x  32  detector 
array  in  (b).  In  both  cases,  the  width  parameters  were  chosen  to  simulate  unresolved  point 
sources. 

The  object  Fourier  domain  statistics  are  straightforward  to  derive  for  the  Gaussian 
binary  star  pairs.  The  mean  object  spectrum  is  [14] 


0(u,v)  =  W‘^amc{Wu,Wv)^hpWpexp{-TT{Wp{u‘^  +  v^))} 
-1-  hgWg  exp{-7r{ws{u^  +  u^))}]  , 


(5.2) 
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(a)  (b) 


Figure  5.1  Sample  binary  star  object  realizations  {hp  =  1,  hs  =  0.5  pixels),  (a)  16  x  16 
detector  (Wp  =  Ws  =  0.5  pixels),  (b)  32  x  32  detector  {Wp  =  Wa  =  I  pixel). 

where  the  sine  function  [22]  is  the  Fourier  transform  of  the  uniform  probability  density 
function  associated  with  a  rectangular  support  constraint  of  dimension  W.  The  correlation 
of  (it,  v)  and  (it',  v')  object  Fourier  components  can  be  expressed  as  [14] 


Rooiu,v]u',v')  =  W‘^sm.c{W{u  —  u')^W{v  —  v')) 

X  [hpiOp  exp{-7rit;p(it^  +  -\-v^  +  v''^)} 

+  exp{— 7rit;5(w^  +  +  v^  +  v'^)}] 

+  W'^sinciWu,  lFi;)sinc*(TFu',  Wv')hphsWpWs 

X  exp{—Tr(wp{u^  +  v^)  +  +  'y'^))} 

+  exp{-7r(«;s(it^  +  v^)  +  Wp{u''^  +  v'^))}]  . 


(5.3) 


Equations  (5.2)  and  (5.3)  provide  the  object  Fourier  domain  statistics  for  all  experiments. 
Each  experiment  trial  begins  by  generating  a  new  random  object  realization  which  is  a 
member  of  the  ensemble  with  these  statistics. 


5.2.2  Measurement  Noise.  Equation  (4.35)  suggests  a  measurement  noise  Fourier 
domain  correlation  matrix  which  is  the  sum  of  two  components.  The  shot  noise  component  is 
the  product  of  the  mean  OTF  and  the  mean  object  spectrum  at  a  difference  frequency  with 
a  scaling  factor  SNR^^.  The  component  associated  with  the  uncorrelated  detector  read  noise 
increases  the  diagonal  elements  by  a  scaling  factor  SNR“^.  Thus,  the  measurement  noise  is 
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correlated  with  respect  to  spatial  frequency.  Since  the  scalar  Wiener  filter  does  not  exploit 
this  information,  the  vector  Wiener  filter  should  have  an  advantage  in  noise  suppression. 

Each  experiment  trial  generates  a  detected  image  which  is  corrupted  by  both  shot  noise 
and  detector  read  noise.  Both  noise  sources  can  be  varied  independently  using  the  K  and 
(7r  parameters  which  correspond  directly  to  light  level  and  detector  read  noise,  respectively. 
An  individual  trial  begins  by  generating  a  normalized  object  realization  and  scaling  to  the 
selected  light  level  as  given  by  K.  A  noiseless  image  is  then  created  by  multiplying  the 
object  spectrum  and  OTF.  A  Poisson  random  number  generator  is  used  to  corrupt  the 
noiseless  image  by  using  its  irradiance  values  as  mean  parameters  in  the  Poisson  distribution. 
Finally,  zero-mean  Gaussian  random  numbers  with  standard  deviation  cTr  are  added  to  model 
signal-independent  detector  read  noise.  This  process  is  repeated  many  times  to  constitute  a 
complete  experiment. 

5.2.3  Performance  Metrics.  To  compare  filter  performance  properly,  metrics  must 
be  used  that  provide  a  realistic  performance  measure.  MSE,  a  correlation  coeflBicient,  and 
MSPE  are  used  to  compare  vector  and  scalar  Wiener  filter  performance  in  this  chapter.  Each 
of  these  performance  metrics  provides  a  dijfferent  method  of  comparison  between  the  true 
and  estimated  object  spectra. 

5. 2. 3.1  Mean  Square  Error.  Consider  the  definition  of  the  error  correlation 
matrix  presented  in  Eqs.  (4.39)-(4.42).  The  diagonal  elements  of  represent  the 
theoretical  MSE  of  the  respective  filter.  This  data  is  presented  in  two  forms  in  this  chapter. 
First,  consider  the  average  MSE  defined  as 

(5.4) 

where  P  is  the  number  of  pixels  in  the  detector  array  as  before.  Second,  consider  a  MSE 
metric  given  as  a  function  of  radial  spatial  frequency  p  via  radial  averaging.  This  quantity 
is  denoted  e^{p).  Radial  averaging  is  the  process  of  averaging  a  two  dimensional  function 
along  concentric  circles  to  produce  a  one  dimensional  plot.  All  MSE  data  in  this  chapter  is 
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normalized  by  the  error  associated  with  the  detected  image.  This  normalization  forces  the 
MSE  data  values  to  be  between  0  and  1  when  the  filters  provide  improved  performance  with 
respect  to  the  detected  image.  Here,  the  subscripts  V,  S,  and  D  denote  MSE  associated 
with  the  vector  Wiener  filter,  scalar  Wiener  filter,  and  detected  image,  respectively. 

The  previous  paragraph  presented  theoretical  MSE  metrics.  Experimental  MSE  data 
can  be  collected  via  Monte  Carlo  simulation.  Here,  the  average  squared  error  per  realization 
is  computed  and  then  accumulated  over  a  large  ensemble  of  images  created  under  identical 
statistical  conditions.  The  average  squared  error  as  calculated  for  the  experiment  trial  is 

4  =  p(0(  -  6i)''(Oi  -  6i).  (5.5) 


The  pertinent  sample  statistics  are  the  mean  of  Eq.  (5.5)  and  variance  of  the  mean  defined 
as  [84] 


1 

- 1) 


1=1 


(5.6) 


where  L  is  the  number  of  images  in  the  ensemble.  The  theoretical  average  MSE  metric 
represented  by  Eq.  (5.4)  and  the  experimental  average  MSE  of  Eq.  (5.5)  were  within 
in  all  cases.  Thus,  the  Monte  Carlo  simulation  used  to  generate  other  performance  metrics 
was  shown  to  perform  accurately  with  respect  to  MSE  [14]. 


5. 2.3. 2  Magnitude  and  Phase  Error.  In  some  cases,  MSE  may  be  a  misleading 
quality  indicator  with  respect  to  the  human  observer  [82].  Clearly,  filter  estimates  can  be 
compared  visually.  However,  visual  comparison  is  only  valid  for  individual  realizations. 
Additional  metrics  are  needed  which  reflect  performance  over  a  complete  ensemble.  In  the 
Fourier  domain,  magnitude  and  phase  play  very  different  roles.  In  fact,  phase  plays  the  more 
important  role  in  many  situations  [61].  The  scalar  Wiener  filter  cannot  compensate  for  phase 
distortions  due  to  noise  which  leads  to  less  deblurring  as  noise  increases  [41].  In  contrast, 
the  vector  Wiener  filter  incorporates  some  a  priori  knowledge  of  the  true  object  phase,  which 
can  help  compensate  for  noise  without  excessive  smoothing.  Thus,  a  separate  comparison  of 
magnitude  and  phase  error  performance  is  needed  to  understand  the  value  of  additional  a 
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priori  object  information.  Also,  the  metrics  given  below  will  be  important  in  demonstrating 
filter  performance  beyond  the  OTF  cutoff  frequency. 

A  correlation  coefficient  between  true  and  estimated  Fourier  spectra  has  been  used 
previously  to  compare  estimation  performance  [92].  In  this  chapter,  the  following  correlation 
coefficient  between  0{u,  v)  and  0(u,  v)  is  used  and  can  be  written  as 

(5.7) 

\/(|0(“.“)I^>(|6(m.»)P) 

|7Q^(M,t;)|  takes  on  values  between  0  and  1.  A  value  close  to  one  implies  the  filter  is 
doing  a  good  job  of  reconstructing  that  Fourier  component.  Now  consider  the  phase  of  the 
reconstructed  object  spectrum.  MSPE  is  defined  as 

v)  =  {{<l>oiu, v)  -  (j)oiu,  v)Y),  (5.8) 

where  (f)o{u,  v)  is  the  phase  of  the  true  object  spectrum  and  (j>d{u,  v)  is  the  phase  of  the  esti¬ 
mated  object  spectrum.  Both  Jqq  and  are  two  dimensional  functions.  Direct  comparison 
of  two  dimensional  functions  is  difficult.  However,  Fourier  spectra  exhibit  a  high  degree  of 
radial  symmetry.  This  symmetry  allows  for  radial  averaging  of  '^qq  and  ip^.  The  data  in 
Sections  5.3  and  5.4  associated  with  these  two  metrics  will  be  presented  in  one  dimensional 
form  using  radial  averaging. 


5.3  Deterministic  OTF 


The  vector  Wiener  filter  derivation  given  in  Eq.  (4.28)  assumed  a  random  OTF  and 
incorporated  a  statistical  description  of  this  quantity.  However,  much  valuable  information 
can  be  gained  from  scenarios  in  which  the  OTF  is  deterministic  and  known  as  shown  in 
Eq.  (4.29).  All  results  presented  in  this  section  incorporate  a  known  OTF.  In  all  cases,  a 
square  pupil  function  is  used  to  compute  the  fixed  OTF  via  Eq.  (2.17).  The  point  spread 
function  (PSF)  associated  with  the  square  pupil  is  a  sinc^(a;,^)  pattern  with  the  width  of 
the  first  zero  crossing  given  by  [22] 


2Xz 


(5.9) 
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Table  5.1  Measurement  Noise  Cases 


K 

SNRk 

SNRr 

10000 

100 

278 

1000 

32 

28 

750 

27 

21 

500 

22 

14 

250 

16 

7 

100 

10 

3 

where  Al  is  the  width  of  the  first  zero  crossing,  A  is  the  imaging  wavelength,  z  is  the 
observation  distance,  and  D  is  the  square  pupil  dimension.  To  study  the  effect  of  the  OTF 
on  filter  performance,  the  OTF  cutoff  frequency  is  changed  by  adjusting  D.  Since  it  is  more 
convenient  to  express  these  quantities  in  a  normalized  pixel  space,  Eq.  (5.9)  becomes 

(5.10) 

where  Alp  is  the  width  of  the  first  zero  crossing  in  pixels,  N  is  the  length  of  one  side  of  the 
detector  in  pixels,  and  Dp  is  the  square  pupil  dimension  in  pixels.  In  all  cases.  Dp  will  be 
used  to  identify  the  pupil  size  and  corresponding  PSF-OTF  for  each  case. 

Table  5.1  provides  corresponding  K,  SNRfc,  and  SNRr  values  for  cases  of  interest  in 
this  section.  The  detector  size  here  is  16  x  16  pixels  with  detector  read  noise  fixed  at  (Jr  = 
electrons  per  pixel  in  all  cases.  K  =  10000  photoevents  represents  high  light  level  and  low 
noise.  Note  that  SNRfc  and  SNRr  provide  an  indication  of  the  relative  contribution  of  the  two 
noise  effects  in  each  case.  For  instance,  when  K  =  10000  photoevents,  SNR^  is  much  lower 
than  SNRr,  which  implies  detector  read  noise  will  not  have  a  large  impact  on  performance. 
When  K  =  100  photoevents,  SNRr  is  lower,  indicating  detector  read  noise  will  play  a  more 
significant  role.  All  MSE  plots  in  this  section  were  generated  via  Eq.  (5.4).  The  l7oo(p)| 
and  p^(p)  plots  were  generated  via  Monte  Carlo  simulation  where  L  =  100,000  images. 

As  presented  in  Section  5.2,  W  represents  the  square  dimension  in  the  center  of  the 
image  plane  where  the  random  binary  Gaussian  object  components  are  allowed  to  exist. 
Thus,  W  represents  the  amount  of  randomness  associated  with  the  object  ensemble.  Figure 
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5.2  compares  scalar  and  vector  Wiener  filter  normalized  average  MSE  performance  as  a 
function  of  W  when  the  pupil  size  is  Dp  =  6  pixels  and  K  =  10000  photoevents  (SNR^  =  100, 
SNRr  =  278).  First,  note  that  the  vector  filter  provides  lower  MSE  over  a  large  range  of  W 
values.  This  improved  performance  is  due  in  large  part  to  the  a  priori  object  Fourier  domain 
correlation  information  used  by  the  new  filter.  The  scalar  filter  has  access  only  to  the  object 
power  spectral  density  (PSD),  which  is  equivalent  to  the  diagonal  values  of  the  Rqo  matrix. 
Since  the  object  ensemble  is  non-stationary,  Rqo  bas  non-zero  off-diagonal  elements.  The 
scalar  filter  cannot  incorporate  this  valuable  information  and,  therefore,  is  no  longer  optimal 
with  respect  to  MSE  as  evidenced  by  Fig.  5.2. 

It  should  also  be  noted  that  the  MSE  associated  with  the  vector  Wiener  filter  increases 
as  the  support  constraint  dimension  W  increases.  Matson  noted  similar  performance  for 
iterative  algorithms  [54] .  His  analysis  showed  that  applying  support  constraints  can  provide 
both  a  superresolution  effect  and  variance  reduction  in  the  noisy  Fourier  data  [54].  Support 
constraints  provide  variance  reduction  by  maintaining  Fourier  domain  correlations  in  the 
data  which  provide  weighted  interpixel  averaging  [54].  The  weights  are  associated  with  the 
Fourier  transform  of  the  support  function.  Thus,  as  support  increases,  the  Fourier  trans¬ 
form  of  the  support  function  narrows,  providing  less  averaging  and  degraded  performance. 
Similarly,  the  vector  Wiener  filter  also  provides  interpixel  averaging  based  on  enforcing  de¬ 
graded  data  Fourier  domain  correlations.  As  W  increases,  support  size  increases  and  the 
Fourier  transform  of  the  support  function  narrows.  The  off-diagonal  elements  of  Rqo  are 
reduced,  which  results  in  a  filter  transformation  matrix  with  less  off-diagonal  structure.  Less 
off-diagonal  structure  in  the  filter  transformation  matrix  is  analogous  to  less  interpixel  av¬ 
eraging  as  each  Fourier  component  is  estimated.  In  the  limit,  when  no  support  constraint  is 
used,  the  filter  transformation  matrix  is  diagonal  and  no  interpixel  averaging  occurs.  Hence, 
the  MSE  performance  of  the  vector  and  scalar  Wiener  filters  is  the  same  for  W  =  16  pixels 
in  Fig.  5.2. 

Figure  5.3  shows  a  single  detected  image  realization  with  the  associated  filter  outputs 
for  W  =  8  pixels,  pupil  size  Dp  =  6  pixels,  and  K  =  10000  photoevents  (SNE^  =  100, 
SNRr  =  278).  Each  mesh  plot  is  normalized  to  a  peak  value  of  unity  for  visual  comparison. 
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Figure  5.2  Normalized  average  MSE,  Syle\,  and  versus  support  constraint  dimen¬ 

sion  W.  Pupil  size  Dp  =  6  pixels  and  K  =  10000  photoevents  (SNRfc  =  100, 
SNRr  =  278). 

Clearly,  the  vector  filter  produces  a  sharper  output  which  more  closely  resembles  the  true 
object  realization,  as  suggested  by  the  MSE  data.  The  mean  object  in  (b)  and  the  PSF  in 
(c)  are  provided  for  the  convenience  of  the  reader  and  later  reference. 

As  noted  in  the  discussion  associated  with  Fig.  5.2,  the  vector  Wiener  filter  can  provide 
superresolution.  Figures  5.4  and  5.5  illustrate  this  idea  by  comparing  vector  filter  |7oo(p)| 
and  (p^{p)  with  those  associated  with  the  scalar  filter.  In  both  plots,  the  pupil  size  is  Dp  = 
6  pixels,  the  support  constraint  dimension  is  W  =  8  pixels  and  K  =  1000  photoevents 
(SNRjfc  =  32,  SNRr  =  28).  Note  that  the  light  level  has  been  reduced,  compared  to  the 
results  depicted  in  Figs.  5.2  and  5.3.  In  addition,  the  spatial  frequency  unity  represents  the 
OTF  cutoff  frequency.  The  scalar  filter  cannot  provide  superresolution  since  it  sets  spatial 
frequencies  beyond  the  OTF  cutoff  to  zero  [41].  Thus,  \'yo6iP)\  scalar  filter  drops 

sharply  beyond  the  cutoff  frequency.  However,  the  vector  Wiener  filter  is  able  to  maintain 
a  high  correlation  coefficient  at  these  same  frequencies.  Figure  5.5  supports  the  idea  of 
a  superresolution  effect  since  (p^{p)  associated  with  the  vector  filter  is  much  lower  beyond 
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Figure  5.3  Normalized  mesh  plots  showing  the  improved  performance  characteristics  of  the 
vector  filter  on  binary  Gaussian  function  objects,  (a)  True  object  realization, 
object  randomness  parameter  W  =  8  pixels,  (b)  mean  object,  (c)  PSF,  pupil 
size  Dp  =  6  pixels,  (d)  detected  image  K  =  10000  photoevents  (SNRfc  =  100, 
SNRr  =  278),  (e)  scalar  filter  estimate,  (f)  vector  filter  estimate. 
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Figure  5.4  Radially  averaged  \'yoo\  as  a  function  of  normalized  spatial  frequency  p  for  the 
scalar  and  vector  filters.  The  spatial  frequency  at  unity  corresponds  to  the 
OTF  cutoff  frequency  for  a  =  6  pixel  pupil  function.  Support  constraint 
dimension  W  =  8  pixels  and  K  =  1000  photoevents  (SNRfc  =  32,  SNRr  =  28). 

unity  spatial  frequency.  As  noted  above,  the  a  priori  object  knowledge  incorporated  in  the 
vector  Wiener  filter  provides  superresolution,  since  W  is  associated  with  a  spatial  domain 
support  constraint.  This  support  constraint  is  incorporated  directly  in  Rqo  and  O.  Support 
constraints  have  a  long  research  history  as  a  means  to  superresolve  data  in  the  Fourier 
domain  [54,55]. 

Now  reconsider  the  form  of  the  vector  Wiener  filter  given  by  Eq.  (4.28)  and  the  scalar 
Wiener  filter  given  by  Eq.  (4.40).  It  is  obvious  that  neither  filter  incorporates  detected  image 
information  beyond  the  OTF  cutoff  frequency.  At  those  high  spatial  frequencies,  the  vector 
filter  relies  exclusively  on  object  statistical  information.  In  contrast,  the  scalar  filter  is  not 
capable  of  incorporating  this  knowledge.  Instead,  the  frequency  components  beyond  the 
OTF  cutoff  are  set  to  zero.  Thus,  the  vector  filter  should  continue  to  perform  better  than 
the  scalar  filter  as  the  OTF  cutoff  frequency  is  adjusted  lower.  Figure  5.6  shows  filter  MSE 
performance  as  a  function  of  the  pupil  size  Dp  when  W  =  8  pixels  and  IC  =  1000  photoevents 
(SNRfc  =  32,  SNRr  =  28).  Figure  5.6  shows  that  the  vector  filter  provides  lower  MSE  than 
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Normalized  Spatial  Frequency,  p 


Figure  5.5  Radially  averaged  MSPE,  as  a  function  of  normalized  spatial  frequency  p 
for  the  scalar  and  vector  filters.  The  spatial  frequency  at  unity  corresponds  to 
the  OTF  cutoff  frequency  for  aDp  —  6  pixel  pupil  function.  Support  constraint 
dimension  W  =  S  pixels  and  K  =  1000  photoevents  (SNRfc  =  32,  SNR^  =  28). 
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Figure  5.6  Normalized  average  MSE,  evle%  and  versus  pupil  size  Dp.  Support 

constraint  dimension  W  =  %  pixels  and  K  —  1000  photoevents  (SNRfc  =  32, 
SNR^  =  28). 

the  scalar  filter  in  all  cases.  As  in  Fig.  5.2,  this  MSE  performance  is  the  combined  result  of 
both  a  superresolution  effect  and  enforcing  detected  data  Fourier  domain  correlations. 

In  Chapter  II,  the  vector  Wiener  filter  was  derived  based  on  the  image  model  given  in 
Eq.  (2.7)  with  the  idea  of  properly  modeling  all  noise  effects.  Thus,  we  are  clearly  interested 
in  the  performance  of  this  new  vector  filter  as  measurement  noise  becomes  more  dominant. 
Figure  5.7  provides  MSE  performance  as  a  function  of  A”  for  W  =  8  pixels  and  pupil  size 
Dp  =  6  pixels.  Clearly,  the  performance  of  both  filters  improves  as  light  level  increases. 
However,  it  must  also  be  noted  that  the  vector  filter  reduces  MSE  by  a  wider  margin  as  light 
level  increases.  At  K  =  100  photoevents  the  vector  Wiener  filter  provides  an  additional  15% 
decrease  in  MSE  below  the  baseline  established  by  the  scalar  filter.  At  K  =  1000  photoevents 
the  decrease  in  MSE  is  37%.  This  trend  is  expected  since  the  off-diagonal  elements  of  the 
measurement  noise  Fourier  domain  correlation  matrix  are  due  to  the  shot  noise.  As 
light  level  is  reduced,  shot  noise  is  less  dominant  and  becomes  more  diagonal,  which 
minimizes  the  importance  of  the  off-diagonal  elements.  As  noted  before,  it  is  access  to  the 
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Figure  5.7  Normalized  average  MSE,  Byle\,  and  e\/e%^  versus  K.  Support  constraint 
dimension  W  =  8  pixels  and  pupil  size  Dp  =  6  pixels. 

ofF-diagonal  components  of  the  correlation  matrices  which  provides  the  difference  between 
the  vector  and  scalar  Wiener  filters.  Figures  5.8  and  5.9  also  illustrate  the  superior  noise 
suppression  of  the  vector  Wiener  filter  for  which  the  support  constraint  dimension  W  and 
pupil  size  Dp  are  unchanged  from  Fig.  5.7.  The  \'iQQ{p)  \  and  metrics  are  superior  across 
all  noise  cases  with  a  superresolution  effect  still  visible  beyond  the  OTF  cutoff  frequency. 

Figure  5.10  provides  detected  image  and  filter  outputs  when  K  =  500  photoevents 
(SNRfc  =  22,SNRr  =  14).  The  scalar  filter  simply  smoothes  the  noisy  data  and,  therefore,  is 
not  able  to  resolve  the  two  object  components  in  (e).  In  contrast,  the  vector  Wiener  filter 
resolves  the  components  in  (f).  Note  how  the  object  Fourier  domain  statistics  enforce  the 
support  constraint  imposed  by  the  parameter  W.  The  next  section  illustrates  vector  Wiener 
filter  performance  when  the  OTF  is  random  and  associated  with  atmospheric  turbulence. 
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Figure  5.8  Radially  averaged  |7ool  as  a  function  of  normalized  spatial  frequency  p  for  the 
scalar  and  vector  filters.  The  spatial  frequency  at  unity  corresponds  to  the 
OTF  cutoff  frequency  for  a  Dp  =  6  pixel  pupil  function  and  support  constraint 
dimension  =  8  pixels.  The  v  and  s  designators  differentiate  between  vector 
and  scalar  filter  traces. 
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Figure  5.9  Radially  averaged  MSPE,  as  a  function  of  normalized  spatial  frequency  p  for 
the  scalar  and  vector  filters.  The  spatial  frequency  at  unity  corresponds  to  the 
OTF  cutoff  frequency  for  a  Dp  =  6  pixel  pupil  function  and  support  constraint 
dimension  PF  =  8  pixels.  The  v  and  s  designators  differentiate  between  vector 
and  scalar  filter  traces. 


5-16 


bject 


0  0  0  0 


(e)  (f) 

Figure  5.10  Normalized  mesh  plots  showing  the  effect  of  measurement  noise  on  filter  per¬ 
formance.  (a)  True  object  realization,  object  support  constraint  dimension 
W  =  8  pixels,  (b)  mean  object,  (c)  PSF,  pupil  size  Dp  =  6  pixels,  (d)  de¬ 
tected  image,  K  =  500  photoevents  (SNRfc  =  22,  SNRr  =  14),  (e)  scalar  filter 
estimate,  (f)  vector  filter  estimate. 


Figure  5.11  Mean  object  used  to  generate  the  random  OTF  data  (W  =  10  pixels). 

5.4  Random  OTF  Due  to  Atmospheric  Turbulence 

In  this  section,  the  OTF  is  random  and  associated  with  atmospheric  turbulence.  Since 
our  objective  here  is  to  study  the  effect  of  turbulence  on  the  vector  Wiener  filter,  the  object 
support  constraint  is  fixed  at  W  =  10  pixels  and  the  images  are  photon-limited  (ar  =  0 
electrons  per  pixel).  Figure  5.11  gives  the  mean  object  associated  with  this  constraint.  All 
data  in  this  section  is  based  on  a  32  x  32  image  array.  A  priori  knowledge  of  the  imaging 
scenario  is  assumed  in  all  cases.  Therefore,  perfect  knowledge  of  the  OTF  statistics  Tid  and 
Rnn  is  available.  However,  no  knowledge  is  assumed  about  the  individual  OTF  realizations. 

Within  the  Monte  Carlo  simulation,  a  random  OTF  is  generated  via  Eq.  (2.17)  using 
the  von  Karman  statistics  and  a  Fourier  series-based  phase  screen  generator  [89].  The  pupil 
plane  residual  phase  aberration  (j)r{x,y)  is  due  to  tilt-removed  distortions.  Thus,  a  rudimen¬ 
tary,  first  order  AO  system  is  simulated.  The  OTF  statistics  were  collected  by  repeatedly 
creating  independent,  random  OTF  realizations  using  the  phase  screen  generator,  removing 
tilt,  and  then  averaging  over  an  ensemble  of  size  10000.  The  optical  system  pupil  is  square 
with  dimension  D.  Fried’s  parameter  rg  is  used  to  indicate  turbulence  strength  [23].  Four 
turbulence  strength  cases  are  examined  using  D/rg  =  1,2,4,  and  a  diffraction-limited  OTF. 
The  ratio  between  the  turbulence  outer  scale  Lg  and  rg  is  fixed  at  Lg/rg  =  100.  Thus,  the 
largest  turbulent  eddies  are  two  orders  of  magnitude  larger  than  the  seeing  cell  size.  Figure 
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5.12  gives  tilt-removed  OTF  statistics  as  a  function  of  the  radially  averaged  spatial  frequency 
variable  p  for  the  cases  noted  above.  The  normalized  spatial  frequency  p  =  1  is  associated 
with  the  diflFraction-limited  cutoff  of  the  optical  system  as  before. 

Recall  the  vector  and  scalar  Wiener  filter  expressions  given  in  Eqs.  (4.28)  and  (4.40). 
Here,  it  can  be  seen  that  the  mean  OTF  plays  a  key  role  in  selecting  the  amount  of  data  infor¬ 
mation  in  the  final  reconstructions.  As  the  OTF  is  attenuated,  performance  will  be  degraded 
as  the  filters  rely  on  the  object  model  to  recover  spatial  frequencies  [12, 13].  Increasing  OTF 
variance  should  also  degrade  filter  performance.  A  metric  is  needed  to  combine  information 
about  the  OTF  mean  and  variance  into  a  single  quantity.  OTF  SNR  is  defined  as 


SNR-h{u,v)  = 


^Vax[H{u,  ?;)] 


(5.11) 


Figure  5.13  gives  SNR'^  as  a  function  of  the  radially  averaged  spatial  frequency  variable  p. 
Notice  when  D/to  =  4,  SNR-^(p)  <  1  across  a  broad  band  of  spatial  frequencies.  In  contrast, 
SNR7^(p)  >  1  out  to  the  diffraction-limited  cutoff  of  the  optical  system  in  the  other  cases. 
This  observation  will  be  useful  in  predicting  filter  performance  with  regard  to  the  results 
given  below. 

As  noted  earlier,  our  primary  objective  is  to  investigate  the  performance  of  the  vector 
Wiener  filter  with  regard  to  atmospheric  turbulence.  However,  some  performance  trends 
associated  with  light  level  should  be  noted.  Figure  5.14  gives  the  vector  Wiener  filter  nor¬ 
malized  average  MSE  versus  K  for  the  four  turbulence  strength  cases.  Both  vector  and 
scalar  Wiener  filter  data  are  presented  for  comparison.  As  expected,  MSE  increases  as  light 
level  decreases  [14].  In  cases  where  more  light  is  available,  the  randomness  associated  with 
atmospheric  turbulence  is  the  dominant  effect.  Also  note  the  performance  degradation  be¬ 
tween  D/ro  =  2  and  4.  This  trend  is  also  expected  based  on  the  increased  attenuation  of  the 
mean  OTF  and  increased  OTF  variance  as  shown  in  Fig.  5.12.  However,  the  relative  size  of 
the  increase  in  MSE  shown  in  Fig.  5.14  seems  to  separate  D/to  =  4  from  the  other  cases. 
Two  questions  naturally  arise  from  this  observation: 
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(a) 


Figure  5.12  Tilt-removed  OTF  statistics  due  to  atmospheric  turbulence  versus  radially 
averaged  spatial  frequency  p  (von  Karman  turbulence  statistics  Lo/tq  =  100, 
D/vo  =  1,2,  and  4).  (a)  Mean,  (b)  variance.  The  OTF  due  to  diffraction  is 
provided  as  a  reference  to  the  reader  in  (a). 
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Figure  5.13 


Normalized  Radial  Frequency,  p 

Tilt-removed  OTF  SNR,  SNR7^(p),  due  to  atmospheric  turbulence  versus  radi¬ 
ally  averaged  spatial  frequency  p  (von  Karman  turbulence  statistics  To/^o  = 
100,  D/to  =  1,2,  and  4).  Note  that  SNR7.{(p)  <  1  when  p  >  0.3  in  the 
D/ro  =  ^  case. 
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Mean  Photoevents  Per  Image,  K 

Figure  5.14  Normalized  average  MSE,  and  eg/e%,  versus  K  for  the  diffraction- 

limited  OTF  and  D/vg  =  1,2,  and  4.  The  v  and  s  designators  differentiate 
between  vector  and  scalar  Wiener  filter  traces. 

1.  Is  there  a  limit  to  how  much  OTF  randomness  the  vector  Wiener  filter  can  handle 
before  performance  is  severely  degraded? 

2.  Can  SNR^^  be  used  to  predict  vector  Wiener  filter  performance? 

To  address  the  questions  posed  above,  consider  the  normalized  MSE  at  individual  spa¬ 
tial  frequencies.  Figure  5.15  shows  radially  averaged  normalized  MSE  for  the  four  turbulence 
strength  cases  under  high  light  level  conditions  {K  =  10000  photoevents).  In  general,  the 
normalized  MSE  associated  with  both  filters  increases  as  p  increases.  This  trend  is  due  to 
the  attenuation  of  the  OTF  and  increased  noise  effects  at  high  spatial  frequencies.  However, 
note  the  jump  in  MSE  that  occurs  when  p  ph  0.2  to  0.4  in  the  D/vo  =  4  case.  Return¬ 
ing  to  Fig.  5.13,  we  observe  that  SNR7^(p)  falls  below  one  when  p  «  0.3  for  this  case.  In 
contrast,  the  other  three  turbulence  strength  cases  exhibit  a  much  more  gradual  increase 
in  MSE  with  p.  Thus,  the  vector  Wiener  filter  seems  to  perform  differently  with  respect 
to  MSE  when  SNR'^(/9)  <  1.  Also,  note  that  the  scalar  Wiener  filter  seems  to  follow  the 
same  trends  except  that  its  performance  is  more  degraded.  The  vector  Wiener  filter  has  a 
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Ci. 


Normalized  Radial  Frequency,  p 

Figure  5.15  Normalized  MSE,  £v{p)/£d{p)  and  £s{p)I£d{p),  for  the  diffraction-limited 
OTF  and  D/tq  =  1, 2,  and  A  (k  =  10000  photoevents).  The  v  and  s  designa¬ 
tors  differentiate  between  vector  and  scalar  Wiener  filter  traces.  Filter  MSE 
performance  is  severely  degraded  at  radial  frequencies  p  where  SNR7^(p)  <  1 
in  the  D/ro  =  A  case. 

particular  advantage  at  high  spatial  frequencies  due  to  interpixel  averaging  provided  by  the 
object  statistical  model. 

Now  consider  a  change  in  light  level.  Figures  5.16  and  5.17  provide  the  same  data  as 
Fig.  5.15  except  light  level  has  been  drastically  reduced  {K  =  1000  and  500  photoevents, 
respectively).  Thus,  MSE  has  increased  across  most  spatial  frequencies.  However,  the  general 
trends  with  respect  to  SNR'^(/9)  remain.  SNR7^(/9)  =  1  is  a  relevant  performance  indicator 
for  this  object  class,  even  at  low  light  level. 

Let  us  now  examine  the  MSPE  performance  of  the  filters  with  respect  to  atmospheric 
turbulence  and  shot  noise.  Figures  5.18  and  5.19  provide  traces  plotted  versus  radially 
averaged  spatial  frequency  p.  The  vector  Wiener  filter  produces  reconstructions  with  less 
MSPE  due  to  superresolution  and  Fourier  data  variance  reduction  associated  with  the  a  priori 
object  statistics.  As  with  the  MSE  above,  note  the  substantial  jump  in  MSPE  in  Fig.  5.18 
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Figure  5.16  Normalized  MSE,  6y{p) / ejj{p)  and  eg{p) / e\,{p) ^  for  the  diffraction-limited 
OTF  and  D/ro  =  1,2,  and  A{K  —  1000  photoevents).  The  v  and  s  designators 
differentiate  between  vector  and  scalar  Wiener  filter  traces. 


Figure  5.17  Normalized  MSE,  e\{p) / e%{p^dind  e%{p)/e%{p),  for  the  diffraction-limited 
OTF  and  D/vo  =  1,2,  and  A{K  =  500  photoevents).  The  v  and  s  designators 
differentiate  between  vector  and  scalar  Wiener  filter  traces. 
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Figure  5.18  MSPE  data  <p^{p).  Turbulence  strength  traces  {K  =  10000  photoevents). 

The  V  and  s  designators  differentiate  between  vector  and  scalar  Wiener  filter 
traces. 

between  turbulence  strength  cases  D/Vo  —  2  and  4.  In  this  plot,  K  =  10000  photoevents, 
which  represents  a  high  light  level.  This  observation  lends  additional  support  to  the  use 
of  SNR’^(p)  as  a  performance  indicator  for  the  vector  Wiener  filter.  Figure  5.19  shows  the 
same  type  of  data  except  light  level  is  varied  and  D/ro  =  2.  Light  level  has  a  minimal  effect 
on  when  K  =  10000  and  K  —  5000  photoevents.  Once  again,  the  vector  Wiener  filter 
provides  superior  performance  across  all  cases.  However,  the  greatest  performance  gains  are 
associated  with  spatial  frequencies  where  SNR7.f(/9)  >1. 

Figure  5.20  shows  detected  image  and  filter  outputs  when  D/tq  =  4  and  K  =  1000 
photoevents.  The  scalar  filter  simply  smoothes  the  atmospheric  turbulence-degraded  data  in 
(e).  In  contrast,  the  vector  Wiener  filter  produces  a  sharper  and  better  resolved  reconstruc¬ 
tion  in  (f).  These  image  realizations  illustrate  a  superior  reconstruction  even  when  filter 
performance  has  been  degraded  due  to  the  strength  of  the  atmospheric  turbulence. 
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Figure  5.19  MSPE  data  ^^(p).  Light  level  traces  {D/ro  =  2).  The  v  and  s  designators 
differentiate  between  vector  and  scalar  Wiener  filter  traces. 

5.5  Summary 

Binary  star  simulation  results  show  that  the  vector  Wiener  filter  provides  superior 
reconstructions  when  compared  to  the  scalar  Wiener  filter  for  non-stationary  object  ensem¬ 
bles.  Comparisons  were  conducted  while  varying  the  object  support  constraint  dimension 
W ,  the  support  of  the  OTF,  and  the  measurement  noise  level.  In  addition,  a  superreso¬ 
lution  capability  of  the  vector  filter  was  illustrated  by  examining  performance  beyond  the 
OTF  cutoff  frequency.  Vector  Wiener  filter  performance  was  also  studied  for  photon-limited 
images  degraded  by  atmospheric  turbulence.  Random  OTFs  associated  with  tilt-removed 
atmospheric  turbulence  cases  were  generated  using  a  Fourier  series-based  phase  screen  gen¬ 
erator  [89].  This  experiment  was  the  first  application  of  complete  OTF  correlation  statistics 
to  the  reconstruction  of  turbulence-degraded  images  [12, 13].  Comparisons  were  conducted 
while  varying  the  turbulence  strength  and  light  level.  The  vector  Wiener  filter  was  superior 
to  the  scalar  Wiener  filter  with  respect  to  normalized  MSE  and  MSPE  across  all  cases  exam¬ 
ined.  However,  substantial  performance  degradation  was  noted  at  spatial  frequencies  where 
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Random  Object 


Mean  PSF 


Scalar  Filter  Estimate 


Vector  Filter  Estimate 


Figure  5.20  Normalized  mesh  plots  showing  the  improved  performance  characteristics  of 
the  vector  Wiener  filter  on  photon-limited  binary  star  objects  degraded  by 
atmospheric  turbulence,  (a)  True  object  realization,  (b)  mean  PSF  {D/to  = 
4),  (c)  random  PSF,  (d)  detected  image  {K  =  1000  photoevents),  (e)  scalar 
filter  estimate,  (f)  vector  filter  estimate. 
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SNR«(/!))  <  1.  Finally,  individual  filter  output  realizations  were  presented  to  demonstrate 
vector  Wiener  filter  capabilities  graphically. 

The  object  irradiance  distribution  was  assumed  to  be  a  random  process  with  known 
first  and  second  order  Fourier  domain  statistics.  Clearly,  the  vector  Wiener  filter  is  not 
applicable  in  situations  in  which  no  information  is  available  about  the  type  of  objects  to  be 
imaged.  The  next  chapter  investigates  the  performance  and  robustness  of  the  vector  Wiener 
filter  for  generalized  object  and  OTF  statistical  models.  The  objective  is  to  study  conditions 
under  which  the  new  filter  is  most  useful. 
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VI.  Vector  Wiener  Filter  Performance  and  Robustness  Study 

6.1  Introduction 

In  Chapter  V,  vector  Wiener  filter  performance  was  investigated  for  an  important  as¬ 
tronomical  object  class,  the  binary  star  pair.  In  that  case,  the  new  filter  outperformed  the 
scalar  Wiener  filter  over  a  wide  range  of  imaging  conditions.  However,  these  conclusions 
are  limited  since  only  binary  stars  and  optical  transfer  functions  representing  atmospheric 
turbulence  were  examined.  In  this  chapter,  performance  and  filter  robustness  will  be  inves¬ 
tigated  for  generalized  object  and  optical  transfer  function  (OTF)  models.  The  objective 
is  to  draw  conclusions  about  the  application  of  the  vector  Wiener  filter  in  general  imaging 
scenarios.  The  following  questions  are  of  primary  interest; 

1.  How  do  the  object  statistics  limit  filter  performance? 

2.  How  do  the  OTF  statistics  limit  filter  performance? 

3.  How  robust  is  the  vector  Wiener  filter  to  errors  in  the  object  statistical  model? 

The  statistical  models  and  resultant  mean  square  error  (MSB)  data  are  based  on  a  16  x  16 
pixel  image  array  in  all  cases. 

This  chapter  is  organized  as  follows.  Section  6.2  presents  the  object  and  OTF  general¬ 
ized  statistical  models.  Section  6.3  reviews  the  filter  expressions  and  relevant  MSB  metrics. 
Section  6.4  gives  performance  data  as  key  object  and  OTF  model  parameters  are  varied. 
Here,  all  statistical  models  are  assumed  to  be  accurate  and  known.  Section  6.5  gives  ro¬ 
bustness  data  when  error  is  introduced  to  key  object  model  parameters.  Finally,  Section  6.5 
summarizes  the  chapter  and  highlights  the  major  conclusions. 

6.2  Generalized  Models 

To  pursue  this  study,  object  and  OTF  statistical  models  are  used  which  can  represent 
a  variety  of  imaging  scenarios.  These  models  are  based  on  parameters  that  control  the  mean 
and  covariance  of  the  respective  random  process.  The  discussion  below  defines  object  and 
OTF  generalized  statistical  models,  key  parameters,  and  supporting  assumptions. 
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6.2.1  Object.  An  object  statistical  model  for  use  in  the  vector  Wiener  filter  is 
defined  by  the  Fourier  domain  mean  and  covariance,  denoted  0{u,v)  and  Coo{UjU',u',v'), 
respectively.  However,  Fourier  domain  quantities  can  be  difficult  to  visualize.  Thus,  the 
generalized  object  model  used  in  this  study  is  defined  via  the  image  domain  statistics  o{x,  y) 
and  Cooix,y,x',y').  These  quantities  are  then  transformed  to  the  Fourier  domain  using  the 
vector-matrix  expressions  [64] 

0  =  Fo,  (6.1) 


and 

Coo  =  FCooF^,  (6.2) 

where  O  and  o  are  P  x  1  ordered  column  vectors  representing  Fourier  and  image  domain  mean 
objects,  respectively;  Coo  and  Coo  are  P  x  P  Fourier  and  image  domain  object  covariance 
matrices,  respectively;  F  is  the  P  x  P  Fourier  transformation  matrix;  H  denotes  a  matrix 
Hermitian  transpose;  and  P  is  the  total  number  of  pixels  in  the  image  array.  An  arbitrary 
element  of  the  matrix  F  is  defined  as  [64] 


Fix,  y\  u,  u)  =  exp  |-^(ux  +  U2/)| ,  (6.3) 

where  a,n  N  x  N  square  detector  array  is  assumed  and  N  =  y/P. 

With  the  relationship  between  image  and  Fourier  domain  statistics  defined,  the  follow¬ 
ing  important  question  remains:  “What  is  a  representative  object  model?”  In  general,  the 
mean  object  provides  low  spatial  frequency  information  about  the  object  class.  A  Gaussian 
function,  centered  at  the  origin,  can  represent  this  low  pass  characteristic.  As  noted  in  the 
previous  chapter,  support  constraints  have  been  widely  used  in  image  reconstruction  algo¬ 
rithms  [41,54].  A  mean  object  based  on  a  Gaussian  function  and  incorporating  a  support 
constraint  can  be  written  as 


oix,y) 


exp  {-■ TT  ((x^  -I-  y^)/alwl)}  (x,  y)  e  support  region 

0  else 


(6.4) 
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(a)  (b) 

Figure  6.1  Mesh  plots  showing  mean  object,  o{x,y),  for  mean  object  width  parameter, 
Wo  =  0.5.  Support  constraint  dimension  (a)  W  =  10  pixels,  (b)  W  =  16  pixels. 

where  the  support  region  is  aW  xW  pixel  rectangular  area  in  the  center  of  the  image  array, 
do  is  a  mean  object  normalization  constant,  and  Wg  is  the  mean  object  width  parameter. 
The  mean  object  normalization  constant  was  chosen  such  that  the  e“^  width  of  the  Gaussian 
function  is  N/2  pixels  when  Wq  =  0.5.  In  addition,  the  mean  object  is  normalized  such  that 


X^o(a;,y)  =  1000,  (6.5) 

x,y 

which  is  associated  with  an  average  photon  count  K  =  1000  photoevents  and  a  moderate 
light  level  as  noted  in  Chapter  V.  Figure  6.1  illustrates  this  mean  object  model  for  Wo  =  0.5. 
The  support  constraint  dimension  is  W  =  10  pixels  in  (a)  and  W  =  16  pixels  in  (b). 

Clearly,  the  object  covariance  will  vary  widely  depending  on  the  particular  imaging 
problem.  However,  standard  image  covariance  models  do  exist  [41].  The  object  covariance 
used  in  this  study  is  based  on  a  two  dimensional,  first  order  Markov  model  and  is  defined  as 


Coo{x,y\  x',y') 


(^o(x,  y)cTo(x',  y')p^  {x,  y)  and  (x',  y')  e  support  region 

0  else 


(6.6) 


where  ao{x,  y)  is  the  object  standard  deviation  at  the  pixel  (x,  y),  po  is  the  object  correlation 
coefficient,  and  d  =  J{x  —  x'Y  +  (y  —  y'Y  is  the  Euclidean  distance  between  arbitrary  pixels. 
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Table  6.1  Object  Model  Parameters 


Parameter 

Description 

Default 

High 

Low 

W 

Support  Dimension 

■■ 

16 

N/A 

Wo 

Mean  Width 

0.9 

0.1 

SNRo 

Image  Domain  SNR 

10 

0.5  or  2 

Po 

Correlation  Coefiicient 

0.5 

0.9 

0.1 

The  object  standard  deviation  ao(x,y)  is  determined  based  on  an  image  domain  signal-to- 
noise  ratio  (SNR)  metric  such  that 


(^oix,y) 


o(x,y)/SNRo  {x,y)  E  support  region 
0  else 


(6.7) 


where  SNRo  is  assumed  to  be  uniform  across  all  object  pixels.  The  incorporation  of  a  support 
constraint  in  the  model  ensures  that  the  object  covariance  is  non-stationary.  The  object 
covariance  is  also  controlled  by  the  SNR^  and  Po  parameters.  As  SNRq  increases,  the  model 
provides  more  precise  information  about  the  object  realizations.  As  po  — ^  1,  the  object  pixels 
are  more  highly  correlated.  Table  6.1  lists  the  key  object  model  parameters  including  values 
which  represent  high,  low,  and  default  cases  in  the  performance  and  robustness  studies. 


6.2.2  Optical  Transfer  Function.  The  object  statistics  described  in  the  previous 
section  are  directly  dependent  on  model  parameters  as  listed  in  Table  6.1.  Once  the  object 
model  parameters  are  selected,  the  mean  and  covariance  are  transformed  to  the  Fourier 
domain  as  described  in  Eqs.  (6.1)  and  (6.2).  A  similar  model  is  also  required  for  the  OTF. 
In  this  case,  the  statistics  will  be  defined  directly  in  the  Fourier  domain.  Both  the  mean  and 
covariance  elements  are  assumed  to  be  real  numbers.  This  assumption  is  consistent  with  an 
OTF  ensemble  representing  imaging  through  atmospheric  turbulence  [74] . 

The  mean  OTF  is  a  low  pass  filter  which  attenuates  high  spatial  frequencies,  yielding 
a  blurred  image.  As  noted  above,  this  characteristic  can  be  represented  using  a  Gaussian 
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Figure  6.2  Mesh  plot  showing  the  mean  OTF,  l-Liu^v),  for  mean  OTF  width  parameter, 
w-H  =  0.5. 

function.  Thus,  the  mean  OTF  is  defined  as 

H{u,  v)  =  exp  {  -TT  +  v^)/a^w^')  }  ,  (6.8) 

where  a'H  is  a  mean  OTF  normalization  constant  and  is  the  mean  OTF  width  parameter. 
As  with  ao  above,  was  chosen  such  that  the  e“^  width  of  the  Gaussian  function  is  N/2 
pixels  when  wy,  =  0.5.  Figure  6.2  illustrates  the  mean  OTF  when  Wn  =  0.5. 

The  OTF  covariance  model  is  the  same  as  the  object  covariance  model  without  a 
support  constraint.  Thus,  u')  is  defined  as 

Cnn{u,  v]  u',  v')  =  an{u,  v)an{u',  v')  p^,  (6.9) 

where  o‘'h{u,  v)  is  the  standard  deviation  of  the  OTF  at  the  spatial  frequency  (u,  v),  py,  is  the 
OTF  correlation  coefficient,  and  d  =  ^{u  —  u'Y  +  (u  —  v'Y  measures  the  relative  distance 
between  arbitrary  spatial  frequencies.  In  the  previous  chapter,  the  OTF  SNR,  SNR'w,  was 
used  to  predict  performance  limits  associated  with  processing  binary  star  images  collected 
through  atmospheric  turbulence.  One  objective  of  this  chapter  is  to  study  performance 
further  with  respect  to  the  SNRt.^  metric.  Thus,  SNR^  will  be  used  to  define  the  OTF 
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standard  deviation  in  Eq.  (6.9)  via  a  simple  manipulation  of  Eq.  (5.11)  which  yields 


(t-h{u,v) 


\H{u,  u)| 
SNR^(u,u)' 


(6.10) 


As  illustrated  in  Fig.  5.13,  a  realistic  SNRt^  function  is  infinite  at  (u,  v)  =  (0, 0)  which  refiects 
a  non-random  DC  value.  The  function  then  rolls  off  to  a  smaller  value  at  a  higher  spatial 
frequency  depending  on  the  attenuation  of  'H{u,  v)  and  the  relative  value  of  0-7^  (u,  v).  For  the 
purposes  of  this  study,  it  is  assumed  that  u)  <  1  for  all  spatial  frequencies  («,  v).  This 
assumption  is  consistent  with  a  realistic  random  OTF  realization  as  defined  in  Eq.  (2.17). 
Based  on  these  characteristics  and  assumptions,  SNR^^  is  modeled  via  the  function 


SNRn{u,v)  =  -^^+n{u,v),  (6.11) 

where  (3-h  is  an  OTF  SNR  normalization  constant  and  is  an  OTF  SNR  roll-off  parameter. 
The  OTF  SNR  normalization  constant  /?«  was  chosen  such  that  Eq.  (6.11)  is  unity  at  one- 
half  the  maximum  radial  frequency  when  fin  =  0.5.  Thus,  controls  the  spatial  frequency 
at  which  SNRt^  falls  below  unity  as  illustrated  in  Fig.  6.3. 

Table  6.2  lists  the  key  OTF  model  parameters  including  values  which  represent  high, 
low,  and  default  cases  in  the  performance  and  robustness  studies.  The  default  value  Wn  =  0.5 
represents  moderate  attenuation  of  the  object  high  spatial  frequencies  by  the  mean  OTF. 
OTF  covariance  associated  with  the  default  values  fi-u  =  1.0  and  pji  =  0.5  represent  moderate 
randomness  and  correlation.  The  term  moderate  refers  to  blur  imposed  by  imaging  through 
the  turbulent  atmosphere  under  good  seeing  conditions.  This  information,  as  well  as  the 
object  parameters  in  Table  6.1,  are  an  important  reference  to  the  reader  with  respect  to  the 
performance  data  given  in  Section  6.4.  The  next  section  reviews  the  filter  expressions  and 
relevant  MSE  metrics  used  in  Sections  6.4  and  6.5. 


6.3  Filter  Expressions  and  Mean  Square  Error  Metrics 

As  noted  in  the  previous  section,  the  mean  object  photon  count  is  fixed  at  K  =  1000 
photoevents.  The  impact  of  changing  light  level  is  well  understood,  based  on  the  data  and 
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Figure  6.3  OTF  SNR,  SNR>^,  versus  radial  spatial  frequency,  p.  The  OTF  SNR  roll-off 
parameter,  p,'n,  controls  the  spatial  frequency  at  which  the  SNR>^  function  falls 
below  unity. 


Table  6.2  Optical  Transfer  Function  Model  Parameters 


Parameter 

Description 

Default 

High 

Low 

wn 

Mean  Width 

0.9 

0.1 

SNR  Roll-Off 

1.0 

3.0 

0.1 

Ph 

Correlation  Coefficient 

0.5 

0.9 

0.1 
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discussion  in  Chapter  V.  Thus,  light  level  is  not  a  variable  in  these  studies  and  a  general 
form  of  the  vector  Wiener  filter,  defined  as 

6  =  Roo'H^  {Roo  0  Rnn  +  Rnn}  ^  D,  (6.12) 

is  applicable.  Equation  (6.12)  is  equivalent  to  Eq.  (4.28)  since  i2o„o„  =  Roo/{R^Y-  All  image 
realizations  are  assumed  to  be  photon-limited  for  these  studies.  Thus,  <7^  =  0  electrons  per 
pixel  and  an  arbitrary  element  of  the  Fourier  domain  noise  autocorrelation  matrix  is 

Rnn{u,  V]  u\  v')  —  H{u  —  u',v  —  v')0{u  —  u',v  —  v').  (6.13) 

The  analogous  scalar  Wiener  filter  is  defined  via  Eq.  (4.40)  as  before. 

In  Chapter  V,  two  MSE  metrics  were  used  to  compare  the  performance  of  the  vector 
and  scalar  Wiener  filters.  In  this  chapter,  only  average  MSE  is  used.  This  approach  supports 
efficient  presentation  of  large  amounts  of  data.  Thus,  in  this  chapter,  the  term  MSE  refers 
exclusively  to  average  MSE  as  defined  in  Eq.  (5.4).  The  required  theoretical  error  correlation 
matrices  are  generated  using  Eqs.  (4.39),  (4.41),  and  (4.42).  In  addition,  vector  Wiener  filter 
MSE  is  normalized  with  respect  to  the  scalar  Wiener  filter  and  detected  image  MSE.  As  in 
Chapter  V,  the  subscripts  V,  S,  and  D  denote  MSE  associated  with  the  vector  Wiener  filter, 
scalar  Wiener  filter,  and  detected  image,  respectively. 

6.4  Performance  Study 

The  objective  of  this  study  is  to  establish  quantitative  limits  on  vector  Wiener  filter 
performance  as  key  object  and  OTF  model  parameters  are  varied  to  represent  a  variety  of 
imaging  conditions.  Here,  all  statistical  models  are  assumed  accurate  and  known  without 
error  as  the  parameters  change.  In  other  words,  the  filters  have  knowledge  of  the  true  imaging 
scenario  and  are  adjusted  accordingly.  Section  6.4.1  outlines  the  important  assumptions 
associated  with  this  work  while  Section  6.4.2  presents  both  object  and  OTF  model  data. 
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6.4-1  Assumptions.  There  are  two  additional  assumptions  associated  with  the 
object  and  OTF  statistical  models  which  pertain  to  the  performance  study.  First,  the  mean 
object  width  parameter  is  fixed  at  Wg  =  0.5.  This  assumption  is  valid  since  a  change  in  mean 
object  size  in  the  image  domain  can  be  viewed  as  altering  the  optical  system  magnification. 
Second,  the  mean  OTF  width  parameter  is  fixed  at  wn  =  0.5  for  all  performance  study  data. 
The  role  of  the  mean  OTF  in  filter  performance  is  relatively  straightforward  as  established 
in  Chapter  V.  Thus,  the  second  order  statistical  quantities  are  of  primary  concern  in  this 
study. 


6-4-2  Data.  Based  on  the  assumptions  noted  above,  this  study  will  concentrate  on 
object  and  OTF  statistical  model  parameters  which  impact  the  covariance.  The  important 
object  parameters  are  the  support  constraint  dimension,  W,  the  SNR  parameter,  SNRg,  and 
the  correlation  coefficient,  po-  The  following  plots  will  illustrate  typical  vector  Wiener  filter 
performance  for  a  variety  of  imaging  conditions. 

Figure  6.4  shows  vector  Wiener  filter  MSE,  normalized  with  respect  to  the  scalar 
Wiener  filter  MSE,  as  a  function  of  the  support  constraint  dimension,  W.  Here,  pg  =  0.5 
and  the  OTF  is  non-random  (SNR-^  =  oo  for  all  spatial  frequencies  (u,v)).  The  data  shows 
that  vector  filter  normalized  MSE  increases  as  W  increases  and  SNRp  decreases.  In  Chapter 
V,  it  was  shown  that  a  rectangular  support  constraint  has  an  important  impact  on  vector 
filter  performance  for  binary  star  pairs.  Support  constraints  provide  variance  reduction 
by  maintaining  Fourier  domain  correlations  in  the  data  which  provide  weighted  interpixel 
averaging  [54].  This  data  clearly  illustrates  the  value  of  interpixel  averaging  in  improving 
vector  filter  MSE  when  object  SNR  is  low.  In  the  SNRp  =  2  case,  performance  improvement 
over  the  scalar  filter  is  less  than  10%  when  no  support  constraint  is  used.  When  W  =  8 
pixels,  this  factor  improves  to  more  than  20%.  Also  note  the  importance  of  SNRq  on  vector 
filter  MSE.  Figure  6.4  shows  that  the  vector  filter  provides  over  40%  less  MSE  than  the 
scalar  filter  when  SNRo  =  10  and  no  support  constraint  is  used. 

Figure  6.4  illustrated  the  importance  of  object  support  on  filter  MSE  performance, 
especially  when  W  is  small.  Many  objects  cannot  be  confined  to  a  small  support  region  and 
are  better  represented  by  a  larger  W  value.  Thus,  the  support  constraint  dimension  is  fixed 


6-9 


Figure  6.4  Normalized  average  MSE,  versus  support  constraint  dimension,  W,  Po  = 

0.5.  The  OTF  is  non-random  (SNR’h  =  oo  for  all  spatial  frequencies  {u,v)). 

at  W  =  10  pixels  for  the  remaining  data  plots  to  represent  a  typical  object  class.  Figure 
6.5  gives  normalized  MSE  as  a  function  of  the  object  SNR  parameter,  SNRq.  Once  again, 
the  OTF  is  non-random  (SNR^^  =  oo  for  all  spatial  frequencies  (u,  v)).  First,  consider  the 
MSE  normalized  with  respect  to  the  scalar  filter  in  (a).  Here,  vector  filter  MSE  approaches 
to  within  10%  of  the  scalar  filter  when  po  =  0.1.  Note  that  the  normalized  MSE  begins  to 
decrease  slightly  when  SNRq  <  2  and  po  =  0.1  or  0.5.  When  SNRp  >  2,  both  filters  perform 
well,  with  the  vector  filter  providing  greater  MSE  improvement  as  SNRo  increases.  When 
SNRo  <  2,  MSE  increases  significantly  for  both  filters,  but  the  scalar  filter  performance 
degrades  faster  than  the  vector  filter.  Thus,  the  object  spatial  SNR  provides  a  natural  per¬ 
formance  threshold  in  these  cases.  The  data  in  (a)  also  shows  that  a  high  object  correlation 
coefficient  can  have  a  substantial  effect  on  MSE.  When  po  =  0.9,  the  vector  filter  offers 
approximately  20%  more  MSE  improvement  over  the  scalar  filter  than  in  the  po  =  0.1  or 
0.5  cases.  In  (b),  the  MSE  is  normalized  with  respect  to  the  detected  image.  In  almost  all 
cases,  the  vector  filter  provides  at  least  20%  improvement.  As  in  (a),  the  MSE  decreases 
when  SNRo  <  2  and  po  =  0.1  or  0.5.  Here,  the  detected  image  MSE  is  increasing  much 
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Figure  6.5  Normalized  average  MSE,  (a)  Syfeg,  and  (b)  Sy/e^j,  versus  the  object  SNR 
parameter,  SNRq.  The  object  support  constraint  dimension  is  fixed  at  W  =  10 
pixels.  The  OTF  is  non-random  (SNR^^  =  oo  for  all  spatial  frequencies  {u,v)). 
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faster  than  the  vector  filter  MSE.  In  general,  the  vector  Wiener  filter  provides  the  best  MSE 
performance  when  SNRo  >  2  and  po  is  high.  This  MSE  improvement  can  be  as  much  as  80% 
over  the  scalar  filter  and  90%  over  the  detected  image  when  SNRq  =  10. 

Now  consider  Fig.  6.6  in  which  normalized  MSE  is  presented  specifically  as  a  function 
of  the  object  correlation  coefficient,  po.  As  before,  the  OTF  is  non-random  (SNR^^^  =  oo 
for  all  spatial  frequencies  {u,v)).  The  data  in  (a)  is  normalized  with  respect  to  the  scalar 
filter.  Note  that  po  =  0.9  provides  20-30%  improvement  over  po  =  0.1  regardless  of  the  SNRq 
case.  Also,  the  object  correlation  coefficient  seems  to  provide  the  greatest  MSE  performance 
boost  over  the  scalar  filter  when  po  >  0.5.  In  (b),  the  vector  filter  MSE  data  is  normalized 
with  respect  to  the  detected  image.  The  plot  shows  that  po  has  the  most  impact  on  MSE 
performance  when  SNRq  is  low.  High  object  SNR  reduces  the  relative  contribution  of  the 
object  covariance  function,  Coo(x,  y,  x',  y'),  to  the  object  correlation  function,  Roo{x,  y;  x',  y'). 
Hence,  po  has  little  effect  on  vector  filter  MSE  when  SNRo  =  10  in  (b).  In  general,  low 
object  SNR  and  po  >  0.5  are  conditions  under  which  the  object  correlation  coefficient  has 
the  greatest  effect  on  filter  MSE  performance. 

For  the  object  model  parameters  investigated  above,  the  OTF  was  non-random  (SNRt^  = 
oo  for  all  spatial  frequencies  {u,  v)).  Now  consider  a  random  OTF  and  investigate  filter  MSE 
performance  with  respect  to  the  OTF  model  parameters.  The  key  parameters  in  this  case 
are  the  OTF  SNR  roll-off  parameter,  p,-^,  and  the  OTF  correlation  coefficient,  p-^.  Figure 
6.7  presents  MSE  as  a  function  of  the  OTF  SNR  roll-off  parameter,  In  (a),  the  vector 
Wiener  filter  MSE  data  is  normalized  with  respect  to  the  scalar  filter  with  the  object  model 
paraimeters  fixed  at  the  default  values  {W  =  10  pixels,  SNR,,  =  2,  po  =  0.5).  Note  that 
vector  filter  MSE  performance  improves  greatly  when  p-^  decreases  below  0.5.  In  fact,  the 
vector  filter  MSE  drops  by  approximately  50%  in  this  case.  At  first  glance,  the  data  in  (a) 
would  seem  to  indicate  that  low  OTF  SNR  improves  filter  MSE  performance.  In  reality, 
both  vector  and  scalar  filters  experience  an  increase  in  MSE  as  p-^  falls  below  this  threshold. 
However,  the  scalar  filter  MSE  increases  more  dramatically  than  the  vector  filter  MSE  as 
shown  via  the  raw  MSE  data  in  (b).  Here,  only  the  pn  =  0.5  case  is  shown.  The  vector 
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(a) 


(b) 

Figure  6.6  Normalized  average  MSE,  (a)  ey/e|,  and  (b)  e\/e%,  versus  the  object  correla¬ 
tion  coefficient,  po.  The  object  support  constraint  dimension  is  fixed  at  W  =  10 
pixels.  The  OTF  is  non-random  (SNRt^  =  oo  for  all  spatial  frequencies  {u,v)). 
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Wiener  filter  is  able  to  maintain  much  better  performance  for  low  OTF  SNR  due  to  the  a 
priori  object  information  provided  by  the  statistical  model. 

Now  consider  the  same  type  of  data  in  Fig.  6.8  in  which  the  object  model  parameters 
are  now  less  favorable  to  vector  filter  processing  {W  =  16  pixels,  SNRq  =  0.5,  po  =  0.5)  and 
Pn  =  0.5.  As  before,  the  vector  filter  MSE  data  is  normalized  with  respect  to  the  scalar  filter 
in  (a).  Note  that  the  MSE  still  decreases  for  low  //•«  values.  However,  the  effect  is  minimal 
compared  to  the  previous  plots.  As  shown  via  the  raw  data  in  (b),  the  vector  filter  MSE 
increases  just  as  fast  as  the  scalar  filter  with  decreasing  In  this  case,  the  object  spatial 
SNR  is  below  the  performance  threshold  established  earlier  and  the  support  constraint  has 
been  removed.  Thus,  the  object  model  cannot  compensate  for  the  increased  randomness 
in  the  OTF.  In  general,  the  vector  Wiener  filter  provides  the  best  MSE  performance  when 
p-H  >  0.5.  This  performance  threshold  corresponds  with  the  OTF  SNR  function  falling 
below  unity  at  normalized  radial  frequency  p  =  0.5.  Thus,  this  data  seems  to  be  in  general 
agreement  with  the  binary  star  data  in  Chapter  V.  For  a  typical  object  class,  vector  Wiener 
filter  MSE  increases  dramatically  when  SNRt.^  falls  below  unity  at  the  mid  and  high  spatial 
frequencies. 

In  Fig.  6.7  (a),  changing  the  OTF  correlation  coefficient  seemed  to  have  little  impact 
on  filter  MSE.  This  performance  trend  is  confirmed  in  Fig.  6.9  in  which  normalized  MSE  is 
given  versus  the  OTF  correlation  coefficient,  pj^.  In  (a)  and  (b),  the  vector  Wiener  filter  MSE 
is  normalized  with  respect  to  the  scalar  filter  and  detected  image,  respectively.  Regardless 
of  the  p-H  value,  the  plot  traces  are  almost  flat.  This  performance  is  relatively  consistent 
across  a  variety  of  object  classes.  Thus,  py,  has  minimal  impact  on  filter  performance. 

6.5  Robustness  Study 

In  some  applications,  the  statistical  models  may  only  approximate  the  true  imaging 
scenario.  This  situation  is  usually  caused  by  a  lack  of  a  priori  knowledge  about  the  problem, 
especially  with  regard  to  the  object  model.  The  objective  of  this  study  is  to  investigate  the 
robustness  of  the  vector  Wiener  filter  to  error  in  key  object  model  parameters.  The  required 
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(a) 


Figure  6.7  (a)  Normalized  average  MSE,  Syleg,  and  (b)  raw  average  MSE  data,  Sy,  versus 

the  OTF  SNR  roll-off  parameter,  Object  model  parameters:  =  10 

pixels,  SNRo  =  2,  po  =  0.5. 


6-15 


(a) 


(b) 

Figure  6.8  (a)  Normalized  average  MSE,  eyje^  and  (b)  raw  average  MSE  data,  versus 

the  OTF  SNR  roll-olf  parameter,  [iji.  Object  model  parameters:  W  =  16 
pixels,  SNRo  =  0.5,  Po  =  0.5. 
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(a) 


(b) 

Figure  6.9  Normalized  average  MSE,  (a)  Sy/es,  and  (b)  ey/e%,  versus  the  OTF  correlation 
coefficient,  p^.  Object  model  parameters:  W  =  10  pixels,  SNRo  =  2,  po  =  0-5. 
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theoretical  error  correlation  matrix  can  be  written  as 


ilee  =  Roo  —  ^wRqd  ~  ^OD^w  +  M^RddM^ ,  (6-14) 

where  Roo,  Rod,  and  Rdd  are  based  on  the  true  model  parameters.  The  filter  transforma¬ 
tion  matrix  is  dependent  on  the  corrupted  model  parameters.  The  rest  of  the  section  is 
outlined  as  follows.  Section  6.5.1  gives  the  important  assumptions  associated  with  this  work 
while  Section  6.5.2  presents  the  data. 

6.5.1  Assumptions.  There  are  several  additional  assumptions  associated  with  the 
object  and  OTF  statistical  models  which  pertain  to  the  robustness  study.  First,  it  is  assumed 
that  the  OTF  model  is  known  with  reasonable  accuracy.  Thus,  only  the  object  model 
parameters  are  studied.  This  assumption  is  valid  since  an  ensemble  of  bright  star  images  is 
often  used  to  obtain  an  accurate  estimate  of  the  atmospheric-optical  system  OTF  for  use  in 
linear  deconvolution  [72,74].  These  same  point  source  images  could  also  be  used  to  estimate 
the  OTF  SNR.  Therefore,  the  mean  OTF  width  parameter  and  OTF  correlation  coefficient 
are  fixed  at  =  0.5  and  =  0.5,  respectively.  Second,  the  object  support  constraint 
dimension  and  the  mean  object  width  parameter  are  fixed  at  W  =  10  pixels  and  =  0.5, 
respectively.  Here,  it  is  assumed  that  these  parameters  are  also  known  with  reasonable 
accuracy.  Thus,  this  robustness  study  will  concentrate  on  the  object  SNR  parameter,  SNRq, 
and  the  object  correlation  coefficient,  Po.  These  parameters  are  associated  with  the  object 
covariance.  It  is  this  part  of  the  object  statistical  model  that  may  not  be  readily  available 
in  some  imaging  applications.  Finally,  percent  error  is  used  as  an  independent  variable  in 
the  study.  This  metric  is  defined  as 

Percent  Error  =  x  100,  (6.15) 

where  is  the  corrupted  model  parameter  and  pt  is  the  true  model  parameter.  In  all 
cases  shown  below,  the  true  parameter  values  are  SNRq  =  2  and  Po  =  0.5.  Percent  error  is 
displayed  from  -100%  to  -i-100%  in  all  cases.  When  the  normalized  MSE  is  unity,  the  vector 
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Wiener  filter  provides  the  same  MSE  as  the  scalar  filter  or  detected  image  for  that  error 
magnitude. 

6.5.2  Data.  First,  consider  error  in  both  the  object  SNR  parameter,  SNRq,  and  the 
object  correlation  coefficient,  po.  Figure  6.10  gives  normalized  MSE  as  a  function  of  percent 
error  in  both  SNRo  and  po.  As  in  the  previous  study,  the  vector  filter  MSE  is  normalized  with 
respect  to  the  scalar  filter  in  (a).  Notice  the  substantial  difference  in  performance  between 
the  p-H  =  0.1  and  //•«  =  1.0  and  3.0  traces.  This  difference  is  associated  with  the  OTF 
SNR  performance  threshold  first  established  in  Fig.  6.7.  Clearly,  vector  filter  performance 
degrades  as  error  is  introduced  into  these  parameters.  In  fact,  the  MSE  performance  ratio 
between  vector  and  scalar  filters  is  greater  than  90%  when  the  error  magnitude  is  greater 
than  60%  in  the  =  1.0  and  3.0  cases.  Now  consider  the  vector  filter  MSE  normalized  with 
respect  to  the  detected  image  in  (b).  The  performance  trends  here  are  similar  to  (a)  except 
the  filter  MSE  is  much  larger  when  the  parameters  are  underestimated  in  the  P'H  =  1.0  and 
3.0  cases. 

The  previous  plots  revealed  the  sensitivity  of  the  vector  Wiener  filter  to  simultaneous 
error  in  the  two  key  parameters  which  control  object  covariance.  Figures  6.11  and  6.12  show 
vector  filter  normalized  MSE  for  error  in  only  one  parameter.  In  both  figures,  the  plot  at  (a) 
is  normalized  with  respect  to  the  scalar  filter  and  the  plot  at  (b)  is  normalized  with  respect 
to  the  detected  image.  In  each  plot,  only  the  p-n  —  1.0  case  is  shown.  Clearly,  the  SNRq 
parameter  generates  the  largest  proportion  of  the  MSE  shown  in  Fig.  6.10.  Thus,  SNRo  is 
more  sensitive  to  error  than  po.  In  fact,  Fig.  6.12  shows  that  error  in  the  po  parameter  never 
degrades  vector  filter  performance  such  that  the  normalized  MSE  is  greater  than  unity. 

6. 6  Summary 

In  this  chapter,  filter  performance  and  robustness  were  investigated  for  generalized 
object  and  OTF  models.  The  objective  was  to  draw  quantitative  conclusions  about  the 
application  of  the  vector  Wiener  filter  in  general  imaging  scenarios.  The  performance  study 
involved  examining  vector  and  scalar  filter  MSE  performance  as  key  statistical  model  param¬ 
eters  were  varied.  In  all  cases,  the  model  parameters  were  assumed  known  without  error. 
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Figure  6.10 


Percent  Error  in  SNR(,  and  Po 

(b) 


Normalized  average  MSE,  (a)  eyle\^  and  (b)  ey/e%,  versus  percent  error 
in  the  object  SNR  parameter,  SNRo,  and  the  object  correlation  coefficient, 
Po-  OTF  model  parameters:  w-^  =  0.5,  =  0.5.  True  parameter  values: 

SNRo  =  ‘2,po  =  0.5. 
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Figure  6.11 


(a) 


Percent  Error  in  SNRp 


(b) 


Normalized  average  MSE,  (a)  and  (b)  versus  percent  error 

in  the  object  SNR  parameter,  SNRj,.  OTF  model  parameters:  W'n  =  0.5, 
pu  =  0.5.  True  parameter  values:  SNRp  =  2,  po  =  0.5. 
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Percent  Error  in  po 

(a) 


Percent  Error  in  po 

(b) 

Figure  6.12  Normalized  average  MSE,  (a)  ey/el,  and  (b)  versus  percent  error  in 

the  object  correlation  coefficient,  Po.  OTF  model  parameters:  w-h  =  0.5, 
p-H  =  0.5.  True  parameter  values:  SNRq  =  2,  Po  =  0.5. 


In  contrast,  the  robustness  study  involved  introducing  error  to  the  filter-assumed  object 
statistical  model. 

The  results  of  the  performance  study  showed  that  vector  Wiener  filter  MSE  perfor¬ 
mance  can  be  limited  by  both  object  and  OTF  statistical  models.  An  object  support  con¬ 
straint  is  extremely  advantageous  to  vector  filter  processing  since  interpixel  averaging  is 
enhanced.  The  object  SNR  parameter,  SNRo,  and  the  object  correlation  coefficient,  po,  pro¬ 
vide  a  fundamental  limit  on  filter  MSE.  When  SNRo  <  2  and  po  <  0.5,  the  vector  filter 
provides  only  marginal  improvement  in  MSE  over  that  of  the  scalar  filter.  However,  Po  >  0.5 
can  help  compensate  for  low  object  SNR  in  many  cases.  The  OTF  SNR  can  also  provide 
a  limit  on  filter  MSE  performance  for  some  object  classes.  In  these  cases,  the  OTF  SNR 
roll-off  parameter,  p-n,  must  be  large  enough  to  boost  the  OTF  SNR  above  unity  at  the  mid 
spatial  frequencies.  For  the  data  shown  here,  the  vector  Wiener  filter  provides  the  best  MSE 
performance  when  >  0.5.  However,  the  vector  Wiener  filter  can  continue  to  perform  well 
below  this  OTF  SNR  threshold  if  the  object  SNR  is  high.  Finally,  it  was  shown  that  the 
OTF  correlation  coefficient,  pn,  has  minimal  impact  on  filter  normalized  MSE. 

The  robustness  study  investigated  the  effect  of  error  in  the  filter-assumed  object  model 
parameters.  These  results  showed  that  the  vector  Wiener  filter  is  less  robust  than  the 
scalar  Wiener  filter  with  respect  to  these  errors.  This  effect  was  anticipated,  since  the 
vector  filter  better  exploits  the  real  world  information  about  the  imaging  scenario  provided 
by  more  detailed  statistical  models.  In  general,  simultaneous  error  in  the  SNRp  and  po 
parameters  resulted  in  marginal  performance  improvement  over  the  scalar  filter  when  the 
error  magnitude  was  greater  than  60%.  The  greatest  impact  on  performance  was  associated 
with  error  in  the  object  spatial  SNR. 
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VII.  Conclusions  and  Recommendations 

7.1  Major  Results 

The  major  results  of  this  research  effort  are  the  following: 

1.  Development  of  a  new  constrained  least  squares  (CLS)  algorithm  for  deconvolution 
from  wavefront  sensing  (DWFS)  processing  of  low  light  images. 

•  Computationally  inexpensive.  In  all  cases  examined,  the  Newton-Raphson  itera¬ 
tion  converged  to  a  solution  in  less  than  10  iterations. 

•  Shot  noise  imposes  the  fundamental  performance  limit. 

2.  Derivation  of  a  new  vector  Wiener  filter  incorporating  the  semi-classical  model  of  pho¬ 
toelectric  light  detection. 

•  Yields  superior  reconstructions  with  respect  to  mean  square  error  (MSE)  and 
mean  square  phase  error  (MSPE)  when  compared  to  the  scalar  Wiener  filter  for 
binary  star  objects. 

•  Provides  superresolution  when  the  object’s  Fourier  domain  statistics  are  known 
for  spatial  frequencies  beyond  the  optical  transfer  function  (OTF)  cutoff. 

3.  Quantitative  results  showing  the  performance  and  limitations  of  the  vector  Wiener 
filter  when  applied  to  binary  star  images  degraded  by  atmospheric  turbulence. 

•  First  application  of  second  order  OTF  statistics  between  different  spatial  frequen¬ 
cies  in  a  Wiener  filter. 

•  Filter  MSE  performance  degraded  for  spatial  frequencies  at  which  the  OTF  signal- 
to-noise  ratio  (SNR)  is  less  than  unity. 

4.  Quantitative  results  showing  performance  limits  on  the  vector  Wiener  filter  associated 
with  generalized  object  and  OTF  models. 

•  The  object  SNR  parameter,  SNRo,  and  the  object  correlation  coefficient,  pg,  pro¬ 
vide  a  fundamental  limit  on  filter  MSE.  When  SNRo  <  2  and  po  <  0.5,  the  vector 
filter  provides  only  marginal  improvement  in  MSE  over  the  scalar  filter. 
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•  For  a  typical  object  class,  the  OTF  SNR  must  be  above  unity  for  normalized 
radial  frequencies  p  >  0.5. 

5.  Quantitative  results  showing  the  robustness  of  the  vector  Wiener  filter  with  respect  to 
object  model  error. 

•  Simultaneous  error  in  the  object  SNR  parameter,  SNRq,  and  the  object  correlation 
coefficient,  po,  resulted  in  marginal  performance  improvement  over  the  scalar  filter 
when  the  error  magnitude  was  greater  than  60%. 

•  The  greatest  impact  on  vector  Wiener  filter  performance  was  associated  with  error 
in  the  object  spatial  SNR. 

7.S  Discussion 

The  previous  chapters  introduced  two  new  linear  reconstruction  techniques  which  com¬ 
plement  existing  linear  filters  and  more  intensive  iterative  optimization  schemes.  The  first, 
CLS  incorporating  wavefront  sensing,  is  practical  for  large  image  arrays  and  easy  to  ap¬ 
ply  when  wavefront  sensor  (WFS)  hardware  is  available.  The  algorithm  is  fundamentally 
limited  by  shot  noise  effects  in  the  phase  estimates.  In  addition,  the  CLS  algorithm  tends 
to  underestimate  the  regularization  constant  for  small  data  ensembles  (<  50  images).  The 
second  technique,  a  new  vector  Wiener  filter,  offers  superior  performance  over  the  existing 
scalar  Wiener  filter  for  non-stationary  image  ensembles.  However,  computational  complex¬ 
ity  severely  limits  the  practical  application  of  this  filter,  since  reconstruction  of  an  AT  x  N" 
array  involves  the  inversion  of  an  x  matrix.  For  a  256  x  256  image  array,  this  means 
inversion  and  storage  of  a  65536  x  65536  matrix!  For  the  vector  Wiener  filter  to  be  widely 
applicable,  methods  must  be  found  to  speed  the  computational  process  and  reduce  memory 
requirements.  The  next  section  offers  some  ideas  for  future  work  related  to  these  limitations. 

7.3  Recommendations  for  Future  Work 

Two  primary  areas  remain  to  be  explored  with  regard  to  these  linear  filter  schemes. 
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7.3.1  Alternate  Constraint  Functions.  In  Chapter  III,  the  Fourier  domain  con¬ 
straint  function  C{u,v)  =  1  was  used  in  the  CLS  algorithm  such  that  the  manual  regular¬ 
ization  constant  e  and  the  inverse  of  the  Lagrange  multiplier  7  were  equivalent.  Alternate 
constraint  functions,  such  as  the  two  dimensional  Laplacian  or  a  support  constraint,  should 
be  investigated.  A  different  constraint  function  may  provide  better  algorithm  performance 
for  smaller  ensemble  sizes. 

7.3.2  Sparse  Matrix  Tools.  In  many  cases,  the  object,  OTF,  and  noise  correla¬ 
tion  arrays  may  be  well  approximated  by  relatively  sparse  matrices.  A  sparse  matrix  is  a 
special  class  of  matrix  that  contains  a  significant  number  of  zero- valued  elements  [53].  This 
important  property  leads  to: 

1.  Reduced  memory  requirements  since  only  the  non-zero  entries  and  their  locations  in 
the  original  matrix  need  be  stored. 

2.  Reduced  computation  time  by  eliminating  operations  on  zero  elements. 

MATLAB  supports  sparse  matrix  computations,  including  a  number  of  iterative  methods  for 
solving  simultaneous  linear  equations  [53] .  These  techniques  could  be  applied  to  the  vector 
Wiener  filter  solution  to  expand  the  practical  computational  limits  of  the  filter. 
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Appendix  A.  Derivation  of  Key  Conditional  Expectations 

This  appendix  presents  the  derivations  of  two  important  conditional  expectations  in  Chapter 
IV.  As  an  aid  to  the  reader,  the  following  additional  information  about  the  semi-classical 
model  of  photoelectric  light  detection  is  provided. 


1.  The  random  variable  K,  the  number  of  detected  photoevents  in  an  image  realization, 
obeys  Poisson  statistics  and  is  described  by  the  probability  density  function  (PDF)  [74] 


PK{K]A)  =  ^(^j  j^X{x,y)dxdy^  exp  1^- j  j^X{x,y)  dx  dy^  (A.l) 

where  A  represents  the  area  associated  with  an  individual  detector  element  and  X{x,y) 
denotes  the  rate  function.  The  rate  function  is  proportional  to  the  noiseless  image 

irradiance  i(a;,  2/). 

2.  The  random  arrival  location  of  the  photoevent  {xn,  yn)  has  a  PDF  related  to  A(a;,  y) 
which  can  be  written  as  [74] 


/  \  _  X{Xn^  yn) 

jf^X{x,y)dxdy- 


(A.2) 


The  mean  number  of  photoevents  occurring  in  the  differential  area  dx  dy  is  A(x,  y)  dx  dy 

[74]. 


A.l  Equation  (4-10) 

In  this  section,  Eq.  (4.10)  is  derived  and  repeated  below 


E, 


K 


v)  exp  {-j2'K{u'xn  +  v'yn)} 


n=l 


=  v')0{u,  v)Ol{u\  v').  (A.3) 


The  reader  should  recognize  that  the  left  side  of  Eq.  (A.3)  is  the  expected  value  of  a  function 
of  the  random  variable  {Xn,yn)-  Recall,  that  the  definition  of  such  an  expectation  is  [62] 


Ex[5(a:)]  =  j  g{a)  p^a)  da. 


(A.4) 
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where  a  is  a  dummy  variable  of  integration.  Rewriting  the  left  side  of  Eq.  (A. 3),  using 
Eq.  (A. 4),  yields 


K 


^Xn,yn\K,'H,0  [•]  =  0{u,v)  /  ^exp{-j27r(u'a;„ + 


.n=l 


^  Pxn,yn  (x„,  yn\K,  n,  O)  dxn  dt/„, 


(A.5) 


where  the  •  notation  represents  the  bracketed  expression  in  Eq.  (A. 3).  Since  both  integration 
and  summation  are  linear  operations,  the  order  of  Eq.  (A.5)  can  be  rearranged  such  that 


^Xn,yn\K,HM  =  0(u,v)f^(f  [p^  n,yni.^ri1  vn\K,n,0) 

n=l 

X  exp  {-j2n{u'xn  +  v'yn)}  dxn  dyn) . 


(A.6) 


The  PDF  associated  with  the  photoevent  arrival  location  was  given  in  Eq.  (A. 2)  above. 
Substituting  this  expression  iov  Px^^y^{xn,yn\K, 1-1,0)  into  Eq.  (A.6)  yields 


'^Xn,yn\K,'H,o[* , 


H  S  S'”) 

f  I  ^\^n,  Vn)  dXfi  dyfp,  \J  J 

X  exp  {-j2Tr{u'Xn  +  v'yn)}  dXn  dyn) . 


(A.7) 


The  integral  in  the  parenthesis  is  the  Fourier  transform  of  the  rate  function.  If  A{u,v) 
denotes  the  Fourier  transform  of  X{x,y),  Eq.  (A.7)  can  be  written  as 


K 


I?  r_i  _  0{u,v)  ^  .j\ 


(A.8) 


where  the  complex  conjugate  of  A{u',v')  is  introduced  based  on  the  positive  sign  of  the 
complex  exponential  kernel.  Recall  that  the  mean  number  of  of  photoevents  occurring  in  the 
differential  area  dx  dy  is  \{x,  y)  dx  dy.  Thus,  the  integration  over  this  quantity  in  Eq.  (A.8) 
is  equal  to  the  average  number  of  photoevents  per  image  K.  Also  note  that  A(m,  v)  can  be 
normalized  such  that 


An{u,v) 


A{u,  v) 


K 


(A.9) 
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Thus,  Eq.  (A. 8)  now  becomes 


Kr.,yn\K,HA*\  =  S  0{u,  v)A*(u',  v').  (A.IO) 

n=l 

Finally,  recall  that  the  rate  function  ^(x,y)  is  proportional  to  the  noiseless  image  i(x,y). 
Thus,  in  the  Fourier  domain  we  can  write  [74] 

A(u,v)  =  ?{(u,v)0(u,v).  (A.ll) 

Making  this  substitution  in  Eq.  (A.IO)  and  noting  that  I2n=i  replaced  with  the 

variable  K  yields  the  final  result  as  stated  in  Chapter  IV 

E..,,.|k,h,oW  =  KH'(u\v')0(u,v)0:(u',v').  (A.12) 

A. 2  Equation  (4-20) 

In  this  section,  Eq.  (4.20)  is  derived  and  repeated  below 

^xn,yn,xm,ym\K,n,o  [•]  =  “  K)On{u,  v)0*^{u\  v')H{u,  v)H*{u\  v') 

+  KH{u  —  u',v  —  v')On{u  —  u',v  —  v'),  (A. 13) 


where  the  •  notation  in  Eq.  (A.  13)  represents  the  bracketed  expression  on  the  left  side  of 
Eq.  (A.  14)  above.  Let  us  begin  with  Eq.  (4.19)  repeated  here  as 


E, 


Xfi ^yn  jXm ^ym 


K 


K  K 


H  {-i27r(ux„  -  u'xra  +  vy^  -  v'ym)} 

Ln=l  m=l 


^Xn,yn\K,'H,0 
"h  ^Xn,yn,Xm,ym\K,'H,0 


I]  exp  {-j27r((M  -  u')xn  +  (?;  -  v')ym)} 

L71=1  . _ _ 

\  K  K 

Y.  exp  {-j27r(ua;„  -  u'x^n  +  vyn  -  v'ym)} 


Ln=l  m=l 


(A.  14) 


n^m 


The  derivation  can  be  divided  into  two  parts:  the  n  =  m  and  n  ^  m  terms. 
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A. 2.1  n  =  m  Term.  First,  rewrite  the  n  =  m  term  using  the  general  expectation 
definition  given  in  Eq.  (A. 4)  such  that 


j  j  exp  {-j2Tr{{u  -  u')Xn  +  (u  -  t^')2/n)}j 


X  Px  n,yn{^nj  Vnl^i  o')  dXn  dj/n 


(A.15) 


As  in  A.l  above,  let  us  rearrange  terms,  substitute  the  rate  function  X{x,  y)  for  the  PDF,  and 
recognize  that  the  resultant  integral  is  the  Fourier  transform  of  A(a;,  y),  such  that  Eq.  (A.15) 


becomes 


A{u  -u',v-  v') 

iEni2/n x  >  ‘ 

n~l  ^ 


(A.16) 


Now  writing  Eq.  (A.16)  in  terms  of  the  normalized  Fourier  domain  rate  function  and  replacing 
the  summation  with  the  variable  K  yields 


^Xn,yn\K,n,On=m  =  ^ Aniu  -  u' ,  V  -  v').  (A.17) 

Finally,  replacing  the  rate  function  with  the  OTF  and  object  spectrum  quantities  gives  the 
n  =  m  term 

^xn,yn\K,n,On^,n  =  KH{u  -  v! ,v  -  v')On{u  -  v! ,v  -  v').  (A.18) 

A. 2. 2  n  ^  m  Term.  As  noted  in  Chapter  II,  one  of  the  key  assumptions  asso¬ 
ciated  with  the  semi-classical  model  is  that  the  number  of  photoevents  occurring  in  non¬ 
overlapping  intervals  are  statistically  independent  [74].  Thus,  we  can  split  the  joint  PDF 
Pxn,yn,xm,ym{^n,yn,Xm,ym)  into  two  marginal  PDFs  such  that 

Pxn,yn,Xm.,ym{Xni  Vm  Vm)  —  Pxn,yn{Xn:  yn)Pxm,ym{Xrmym)-  (A. 19) 

Using  this  product  of  marginal  PDFs  and  writing  the  n^m  term  using  Eq.  (A. 4)  gives 
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^Xn,yn,Xm,ym\K,'H,On:^rn  ~ 


-  k) 

X  j  ex^{-j2'K{uXn  +  vyn)} 

^  PxntVni^nj  Vn)  dXn  dyn) 

X  (//exp  0-2  7T{u'Xm  +  V’ym)} 

^  Pxm,ymi^Tni  Vm)  dxjji  . 


(A.20) 


As  before,  substitute  the  rate  function  X{x,y)  for  the  PDFs  and  note  that  the  integrals  are 
Fourier  transforms  which  yields 


^Xn,yn,Xm,ym\K,'H,Onjtm  ~ 


A(u,  v)A*{u',  v') 

W  ' 


(A.21) 


Finally,  normalizing  and  replacing  the  rate  function  with  the  OTF  and  object  spectrum 
quantities  gives  the  n^m  term 


^xn,vnMm\K,H,On^m  =  “  K)On{u,  v)Ol{u\  v')'H{u,  v)H* {u' ,  v').  (A.22) 

Now  combining  Eqs.  (A. 18)  and  (A.22)  gives  the  Chapter  IV  result 

^xn,yn,xm,ym\K,H,o  [•]  =  -  K)0„{u,  v)0*{u',  v')'H{u,  v)H* {u' ,  v') 

+  KH{u  —  u',v  —  v')On{u  —  u\v  —  v'),  (A. 23) 
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Appendix  B.  Detected  Image  Probability  Density  Function 

In  this  appendix,  it  is  shown  that  Pd|o(D|0)  approaches  a  Gaussian  distribution  for  large 
photoevent  counts  K  based  on  the  Fourier  domain  image  degradation  model  given  in  Eq.  (2.7) 
and  repeated  here  as 

K  p 

D{u,  v)  =  Y,  exp  {-j2-K{uXn  +  vyn)}  +  {-j‘2Tr{uXp  +  vyp)}  .  (B.l) 

n=l  p=l 

In  general,  O  is  random.  Therefore,  the  conditional  probability  density  function  (PDF) 
given  a  specific  realization  of  the  random  object  is  considered.  A  given  realization  for  the 
underlying  optical  transfer  function  (OTF)  H  is  also  assumed. 

Let  us  consider  the  second  term  first  and  rewrite  using  Euler’s  identity  [23]  as 

p 

^  Up  exp  {-j27r{uxp  +  vyp)}  = 

p=i 


Now  recall  that  the  detector  read  noise  Up  was  assumed  to  be  a  zero-mean,  uncorrelated, 
Gaussian  random  variable  with  uniform  variance  af.  Thus,  the  real  and  imaginary  parts  of 
Eq.  (B.2)  are  sums  of  scaled  Gaussian  random  variables  of  the  form 


^  WpRe  [exp  {-j2Tr{uXp  +  vyp)}] 

\p=\ 

+  j  Im  [exp  {-j2n{uxp  +  vyp)}] 

\p=i 
p 

UpCos  [-2'K{uXp  -f  vyp)] 

\p-i 

I  ^ 

+  i  l-2n{uXp  -I-  vyp)]  )  . 

Vp=i 


(B.2) 


and 


Re 


np  exp  {-j2TT{uXp  +  vyp)} 


p=i 


—  kitii  -t-  k2n2  +  ...  +  kpUp, 


(B.3) 


Im 


P 

Tip  exp  {-j2'jT{uXp  -1-  vyp)} 

p=i 


—  llTli  l2n2  T  •••  "I"  IpUp, 


(B.4) 
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where  the  real  constants  kp  and  Ip  are  cosine  and  sine  functions,  respectively.  Written  in 
this  form,  both  the  real  and  imaginary  parts  of  the  second  term  of  Eq.  (B.l)  are  Gaussian 
distributed,  since  scaled  sums  of  statistically  independent  Gaussian  random  variables  are 
also  Gaussian  distributed  [30].  Thus,  the  complex  random  variable 

p 

5^npexp{-j27r(MrEp  +  i;?/p)},  (B.5) 

p=i 


is  Gaussian  distributed. 


Now  consider  the  first  term  of  Eq.  (B.l).  Here,  the  Gaussian  nature  is  less  obvious  since 
the  underlying  PDF  associated  with  the  random  arrival  location  (a:„,  j/„)  is  not  Gaussian  [74]. 
As  before,  let  us  rewrite  in  terms  of  real  and  imaginary  parts  using  Euler’s  identity  such  that 
the  first  term  of  Eq.  (B.l)  becomes 


K 


K 


exp  {-j27r(ua;„  +  U2/„)}  = 

n=l 


53  Re  [exp  {-j2'K{uXn  +  vyn)}\ 


\n=l 


K 


+  i  53  +  Vyn)}] 


\n=l 


K 


=  51  Cos  [-27r(ua;„  +  vyn)] 


\n=l 


+  3  ^53  [-27r(ua;„  +  z;y„)]  j  .  (B.6) 


In  Chapter  II,  it  was  noted  that  non-overlapping  photoevent  arrival  locations  are  statistically 
independent  based  on  the  semi-classical  model  [74].  Thus,  the  real  and  imaginary  parts  in 
Eq.  (B.6)  are  sums  of  independent  random  variables.  In  fact,  these  quantities  are  the  result 
of  large  sums  of  independent  random  variables  since  K  is  on  the  order  of  1,000  to  1,000,000 
for  typical  astronomical  imaging  applications.  We  argue  that  the  first  term  of  Eq.  (B.l)  has 
an  approximate  Gaussian  PDF  based  on  the  Central  Limit  theorem  which  can  be  stated  as 
follows  [23,62]: 

Given  N  independent  random  variables  Xi  with  arbitrary  PDFs  (not  neces¬ 
sarily  the  same) ,  we  form  their  sum 

X  =  Xi  +  X2 -i- Xn.  (B.7) 
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The  random  variable  x  has  mean  ry  =  771+772  +  — +^n  and  variance  =  o-f +  cr|  + 

...  +  «7^.  The  Central  Limit  theorem  states  that  under  certain  general  conditions, 

X  approaches  a  Gaussian  distribution,  with  the  same  mean  and  variance,  as  N 
increases. 

When  the  Xi  are  identically  distributed,  the  general  sufficient  condition  for  application  of 
the  theorem  is  that  the  means  and  variances  of  the  random  variables  must  be  finite  [23]. 
The  sums  of  sines  and  cosines  shown  in  Eq.  (B.6)  fit  the  requirements  of  the  Central  Limit 
theorem  in  that  the  means  and  variances  are  guaranteed  to  be  finite  and  K  is  very  large  for 
typical  applications.  Ref.  [62]  states  that  7i  =  30  is  sufficient  in  most  applications  in  which 
the  random  variables  are  independent  and  identically  distributed.  Thus,  the  first  term  of 
Eq.  (B.l) 

K 

5^exp{-j27r(Ma;„  +  7;7/n)},  (B.8) 

n=l 

approaches  a  Gaussian  distribution  for  this  imaging  application. 

Now  combining  the  conclusions  regarding  the  first  and  second  terms  of  Eq.  (B.l),  it 
can  be  seen  that  D(u,  v)  is  the  sum  of  two  independent  Gaussian  random  variables.  Thus, 
Pd|o(D]0)  is  also  Gaussian,  since  the  sum  of  two  arbitrary  Gaussian  random  variables  is 
Gaussian  distributed  [30]. 
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commission  in  June  of  1989.  Second  Lieutenant  Ford  then  served  as  an  electronics  engineer 
at  Wright  Laboratory’s  Armament  Directorate,  Bglin  Air  Force  Base,  Florida.  His  duties 
included  program  management  responsibility  for  the  Ballistic  Sight  Technology  Improving 
Night /Day  Gunnery  program  which  sought  to  improve  the  first  burst  accuracy  of  small  crew- 
served  weapons.  Ls  December  1994,  Captain  Fbrd  graduated  form  the  Air  Force  Institute  of 
Technology  (AFTT)  with  a  Master  of  Science  Degree  in  Electrical  Engineering.  He  entered 
the  Ph.D.  program  at  AFIT  in  1995. 
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